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Abstract 

We introduce elliptic weights of boxed plane partitions and prove that 
they give rise to a generalization of MacMahon's product formula for the 
number of plane partitions in a box. We then focus on the most gen- 
eral positive degenerations of these weights that are related to orthogonal 
polynomials; they form three two-dimensional families. For distributions 
from these families we prove two types of results. 

First, we construct explicit Markov chains that preserve these distri- 
butions. In particular, this leads to a relatively simple exact sampling 
algorithm. 

Second, we consider a limit when all dimensions of the box grow and 
plane partitions become large, and prove that the local correlations con- 
verge to those of ergodic translation invariant Gibbs measures. For fixed 
proportions of the box, the slopes of the limiting Gibbs measures (that 
can also be viewed as slopes of tangent planes to the hypothetical limit 
shape) are encoded by a single quadratic polynomial. 
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1 Introduction 

The uniform distribution on boxed plane partitions (equivalently, lozenge tilings 
of a hexagon) is one of the most studied models of random surfaces. There are 
four principal types of results regarding this model that have been proved. 

(1) Law of large numbers: Under the global scaling (the bounding box/he- 
xagon is fixed and the mesh is going to zero), the measure concentrates on 
surfaces that are close to a certain deterministic limit shape. The limit shape 
can be obtained as the unique solution to a suitable variational problem. The 
solution is encoded by a second degree polynomial in two variables, see [CKP], 
[CLP], [DMB], [Des], [KO]. 

(2) Locally near any point of the limit shape, the measure on tilings converges 
to a (uniquely defined, see [Sh]) translation-invariant ergodic Gibbs measure on 
lozenge tiling of the plane of a given slope, and the slope coincides with the 
slope of the tangent plane to the limit shape at the chosen point, see [Gor] and 
also [Kcl], [Ke2], [KO], [KOS]. 

(3) The correlation kernel of the random point process of lozenges of one of 
the types is explicitly expressed in terms of classical Hahn orthogonal polyno- 
mials, see [Gor], [Jl], [J2], [JN]. 

(4) A few algorithms, both asymptotic and exact, have been proposed to 
generate the random tilings in question, see [BG], [Kr], [PI], [P2], [Wil], [Wi2]. 

These are complemented by the classical MacMahon product formula for the 
total number of plane partitions in a given box, see, e.g., Section 7.21 in [St]. 

In this paper we study measures on boxed plane partitions that generalize 
the uniform distribution. The weight of a tiling is defined as the product of 
certain simple factors over all lozenges of a fixed type, see Section 2.2 for defi- 
nitions. One special case is the weight q volume j where volume is the volume of 
the corresponding plane partition, and q is an arbitrary positive number. 
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In the most general case we consider, the weight of a lozenge is elliptic. 
Our initial motivation came from the fact that these weights lead to a nice 
generalization of the MacMahon formula mentioned above, see Theorem 10.5 in 
the Appendix. 

For the asymptotic analysis, we look at the top degeneration of the elliptic 
weight that is related to orthogonal polynomials. 

Our asymptotic results amount to proving analogs of (2), (3), and (4) above. 
As for (1), we derive the corresponding variational problem (which differs from 
the one for the uniform case by the presence of an external potential) , and show 
that the hypothetical limit shape (obtained from (2)) solves the corresponding 
Euler-Lagrange equation. However, we do not prove the concentration phe- 
nomenon rigorously. 

One interesting feature of the limit shapes that arise is that the curve that 
bounds the frozen regions may have one or two nodes in vertices of the hexagon, 
see Section 9 for illustrations. 

In terms of orthogonal polynomials, our models go all the way up to the top 
of the g-Askey scheme — the most general models we analyze asymptotically 
are related to the q-Racah classical orthogonal polynomials. We also show that 
the elliptic weights lead to the biorthogonal functions constructed in [SZ]. We 
hope to return to the asymptotic analysis of this case in a later publication. 

Our proof of (2) follows the same steps as the argument for the uniform case 
in [Gor] . It is based on the general method of computing limits of correlation 
kernels suggested in [BO] and [O]. The crucial property we need is the second 
order difference equation satisfied by the q-Racah orthogonal polynomials. 

Our perfect sampling algorithm is obtained from a more general construction 
of relatively simple Markov chains that change the size of the box (one side 
increases by 1 and another side decreases by 1) , and that map the measures from 
the class we consider to similar ones. The construction follows the approach of 
[BF]; the key facts that make that approach possible reduce to certain recurrence 
relations for the q-Racah polynomials. 

A computer simulation of the above-mentioned Markov chains can be found 
at http : //www. math. caltech.edu/papers/Borodin-Gorin-Rains . exe. 

Acknowledgements. AB was partially supported by NSF grant DMS-0707163. 
VG was partially supported by the Moebius Contest Foundation for Young Sci- 
entists. EMR was partially supported by NSF grant DMS-0833464. 

2 Model and results 

2.1 Combinatorial interpretations 

For any integers a, b, c > 1 consider a hexagon with sides a, b, c, a, b, c drawn 
on the regular triangular lattice. Denote by £l a xbxc the set of all tilings of 
this hexagon by rhombi obtained by gluing two of the neighboring elementary 
triangles together (such rhombi are called lozenges). An element of ^3x3x3 is 
shown in Figure 1. 
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Figure 1. Tiling of a 3 x 3 x 3 hexagon. 



Lozenge tilings of a hexagon can be identified with 3-D Young diagrams 
(equivalently, boxed plane partitions) or with stepped surfaces. The bijection 
is best described pictorially. We show a 3-D shape corresponding to a tiling in 
Figure 1. 

It is convenient for us to slightly modify both hexagon and lozenges by means 
of a simple affine transform of the plane. 




Figure 2. Affine modification of lozenges 



We thus obtain a tiling of a slightly different hexagon. 




T=b+c 

Figure 3. Modified tiling of a 3 x 3 x 3 hexagon and the corresponding family 

of non-intersecting paths. 

In what follows we use different parameters instead of a, b, c. Set N = a, 
T=b + c, S = c. We will also denote the set ft axbxc by Q(N, T, S). 
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Each tiling corresponds to a family of nonintersecting paths as shown in 
Figure 3. 

Consider a section of our family of paths by a vertical line t — tg. Clearly 
this gives an iV-tuple of points in Z. Thus, our tiling can be viewed as an 
TV-point configuration varying in time t = 0,1,..., T. Note that when t 
( ) the configuration consists of points {0, 1, . . . , N — 1}, while for t = T the 
configuration consists of points {S, . . . , S + N — 1}. 



2.2 Probability models 

Let us introduce the probability measures on fl(N, T, S) that are studied in this 
paper. For any T e fl(N,T, S), we define its weight w(T) and consider the 
probability distribution on f2(iV, T, S) given by the formula 

Prob{T} = - W{T) 



The weights we consider are such that the probability of a tiling is propor- 
tional to the product of certain weights corresponding to the horizontal lozenges 
in it, i.e., 

W (T)= n «,(<>) 

Note that the number of horizontal lozenges in a tiling of a given hexagon 
is fixed (i.e., it does not depend on the tiling). Hence, multiplying w (' - -') by a 
nonzero constant does not change the probability distribution. 

In the most general case considered in Section 10, w( < C>) is an elliptic weight 
given by 

(<>) = Knal^V-'^p^- 1 !!^) m 

W[ ' ' fl p (gI- a */a-l Ul , gJ-3</2 Ul( gi+3i/2-l U2) gJ+3i/2 U2 ) ' ^ > 

where the coordinates of the topmost point of ' _> are (the i and j axes are 
pictured in Figure 1), u\, 112, p, q are (generally speaking, complex) parameters, 



e p (x) = l[(i- P i x)(i-p i + 1 /x) 



i=0 

and 8 p (a, b,c. . . ) = 9 p {a)6 p {b)9. p {c) . . . 

Mostly we will not be concerned with this most general case, nor with the 
most general trigonometric case obtained by taking p — > (note 6q(x) — 1 — x), 
as in these cases the kernels involve biorthogonal functions. The most general 
orthogonal polynomial case is the limit 

p^O, ui = 0(y/p), u 2 ^0(^/p) 1 uiu 2 = pK 2 q~ s , (2) 
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in which case the weight function is 

« = ^ (S+1)/2 -^ ® 

this is also the most general case in which the weight is independent of i. 

Clearly, the factor q~( s+1 )/ 2 can be removed if we replace n by k! = k ■ 
q~i s + 1 )/ 2 . However, this choice is more convenient for our further considera- 
tions. 

We need to make sure that the probabilities of tilings are positive. This 
leads to certain restrictions on the parameters. There are three possible cases: 

(i) . imaginary q-Racah case: q is a positive real number, k is an arbitrary 

pure imaginary complex number; 

(ii) . real q-Racah case: q is a positive real number, k is a real number with 

additional restrictions depending on the size of a hexagon; n cannot lie 
inside the interval [q~ N+1 ^ 2 , gr( T_1 )/ 2 ] or [q ( - T ~ 1 '>/ 2 , q~ N+1 ^ 2 ], depending 
on whether q > 1 or q < 1; 

(iii) . trigonometric q-Racah case: q and k are complex numbers of modulus 1, 

q = e la , k = , plus additional restrictions on k depending on the size 
of a hexagon: both —a(T— l)/2 + /3 and a(N — 1/2) + [5 must lie in the 
same interval of the form [nk, ir(k +1)], k e Z. In this case 



Kq 



j-(S+l)/2 



2x/=Tsin (a(j-(S + l)/2) + 0), 



Kq j-(S+l)/2 

and the factor 2^/^l here can be omitted. 

The names of the cases are related to those of the classical orthogonal poly- 
nomials that appear in the analysis. 

Denote the resulting measure on fl(N, T, S) by n(N, T, S, q, k). 

There are further limit transitions. 

If we send n — > then we get the q-Hahn case 

w(0) = q~ j . 

Thus, the probability of the plane partition of volume (=numbcr of 1 x 1 x 1 
boxes) V is proportional to q~ v . On the other hand, if we send k — > oo we 
will get the weights q v . Therefore, the case of general k can be viewed as an 
interpolation between the measures q volume and q~ volume . 

Another possibility is to set k = q K and then send q — > 1. The weight of a 
horizontal lozenge tends to 

w(0)= K + j-(S+ l)/2. 

We call this case the Racah case. One has to impose restrictions to ensure 
positivity: K cannot lie inside the interval [— N + 1/2, (T — l)/2]. 
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Finally, if wc cither send k — > and set q = 1, or send g — > 1 and then send 
— > oo, then we obtain the Hahn case 

w(0) = 1 

In this case our probability distribution on f2(iV, T, S 1 ) is uniform. 

Below we mostly work with the imaginary g-Racah case, but all the results 
can be carried over to all the other cases mentioned above by the appropriate 
substitutions of parameters and degenerations. 



2.3 Representation of the weight 

Let us view our tiling as a pile of 1 x 1 x 1 cubes in the box located between the 
planes x\ = 0, x\ = a, x 2 = 0, x 2 = b, x 3 = 0, X3 = c. Then Figure 1 represents 
a projection of the border of the 3-D-diagram to the plane x\ + x 2 + x% = 
parallel to the vector (1, 1, 1). 

For any v £ K 3 , denote by h(v) the Euclidian distance from v to the plane 
X\ + x 2 + X3 = divided by v3, and denote by h(v) the distance from v to the 
union of coordinate planes x\x 2 x^ — computed along the (1, 1, Indirection 
and divided by \/3- 

Recall that the weight of a tiling is given by the formula: 

w{T) = const ■ j~| w(j), 



where j is as in Figure 1. 

Grouping all the lxlxl cubes of the plane partitions into columns with 
fixed coordinates (x 2 , £3), we can rewrite the weight in the form 

w(T) = const ■ J| W{3) 



w(j - 1) 



where the product is taken over all cubes of the plane partition, and j stands 
for the j-coordinate (see Figure 1) of the top vertex of the cube. 
Collecting factors with the same j, we obtain 



h{v)-h(v) 



where the product is taken over all points v on the border of the plane partition, 
whose three coordinates [xi, x 2 , X3) are integers. Equivalcntly, one can think 
of the product being taken over all vertices of the triangular lattice inside the 
hexagon. Note that we may replace h(v) — h(v) by h(v) since the remaining 
product depends only on (N, S, T) . This gives 
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For the g-Racah case we obtain 



w(T) — const ■ 



while for the g-Hahn case 



Kq 



j-(S+l)/2 



1 



h(v) 



Kq 



j-(S+l)/2 



Kq 



j-(S+3)/2 _ 



1 



K gi-(S+3)/2 / 



w(T) = const ■ JJ 



7-J + l 



const • J| g~ h(u) = const ■ q~ v 



oiume 



2.4 Results and variational interpretation 

Our results are of two kinds. 

First, for each of the probability distributions on tilings described above 
we construct explicit discrete time Markov chains that relate random tilings of 
hexagons of various sizes. The elementary steps of these chains change the size 
of the hexagon from a x b x c to a x (6 =p 1) x (c ± 1). 

Randomness in each step consists in generating finitely many independent 
one-dimensional random variables. It takes 0(a(b + c)) arithmetic operations to 
produce a tiling of the ax (6—1) x (c+1) hexagon using a tiling of the axbx c 
hexagon. 

Together with the trivial observation that there is exactly one tiling of an 
a x (6 + c) x hexagon, these chains provide, in particular, an efficient per- 
fect sampling algorithm for random tilings distributed according to the real 
g-Racah, imaginary g-Racah, trigonometric g-Racah, g-Hahn, Racah and Hahn 
distributions. 

A description of the algorithm can be found in Section 6, and in Section 9 
we provide some pictures generated using this algorithm. 

Second, we evaluate the asymptotics of the local behavior of our measures 
as all sides of the hexagon tend to infinity comparably. 

It is known, see [Kel], [Kc2], [OR], [Sh], [KOS], [BS], that for any three posi- 
tive numbers (pi,P2,P3) with P1+P2+P3 — 1, there exists a unique translation- 
invariant Gibbs measure on lozenge tilings of the whole plane such that in a 
large box, the numbers provide asymptotic ratios of the number of 

lozenges of three types. It is convenient to encode the triple (pi,P2,P3) by a 
complex number z in the upper-half plane so that the triangle (0, 1, z) has an- 
gles (TTpi,irp2,iTP3)- The correspondence between angles and lozenge types is 
indicated in the figure below. 
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If one of the pj's tends to 1 (for example, z tends to a point in R away from 
{0,1}), then the corresponding measure degenerates to the "frozen" (nonran- 
dom) tiling of the plane by lozenges of the corresponding type. 

Theorem 2.1. Introduce a small parameter e -C 1, and set 

S=Se- 1 + o(e- 1 ) : T = T£- 1 + o(e- 1 ), TV = Ne -1 + o(e _1 ), q = q £+ °( £ ). 

T/ien near a given point (t,x) = (te _1 ,xe _1 ), the random tiling converges, as 
e — > 0, to a certain ergodic translation-invariant Gibbs measure. The parameter 
z of this measure (encoding the slope (pi,P2iP3) via angles as described) is the 
unique solution in the upper half-plane to the (quadratic) equation 

Q(u,v)=0, (4) 

where 



zq l - K 2 q- s+2x (1 - z)q x 



1 — zn 2 q s+2x * ' 1 — ZK 2 q~ s+2x ' 

and Q is a degree 2 polynomial in (it, v): 



(5) 



Q(u, v) = u 2 

q T - S - N + K 2 (l + q-S+N+T + q -2S+T + q -S-N _ q -S _ q -S+ T) + K 4 q -S+N^ y 2 

\ T - S + q- N + K 2 (q N + q- S )) UV (q T + q 7 ^ + n\l + q^ 1 )^ v 

-(l + q T )u + q T . (6) 

If the solutions to this equation in z are real, then the limit measure is frozen. 



Let us now explain how one could guess these formulas. In Section 8 we 
present a rigorous proof of Theorem 2.1 which uses an argument of a different 
kind. The remainder of this section is purely empirical; we hope to address the 
same issues rigorously in a later publication. 

Although Theorem 2.1 describes the microscopic behavior of our model, the 
parameters (pi,P2,Ps) of the limit measure are closely connected with macro- 
scopic properties. 

If we view tilings as stepped surfaces in a box and scale them in such a way 
that the bounding box remains fixed, then it is plausible that in the limit we will 



9 



observe a deterministic limit shape. The normal vector to this limit shape at 
any point has to coincide with the vector (pi,P2,P3) of the local limit measure 
at this point. 

The concentration of the measure near the limit shape is known in the q- 
Hahn case, see [CKP], [KO]. The limit shape is the unique solution of a certain 
variational problem. It is not hard to pose such a variational problem in the 
g-Racah case as well. 

Recall that in Section 2.3 we found the following representation for the weight 
of a tiling: 



;(T) = const ■ Y\ 



f Kq j-(S+i)/2 _ 



Kq- 



1 

i-(S+l)/2 



\ 



h(v) 



Kq- 



]-(S+3)/2 _ 



K qj-(S+3)/2 ) 

Taking the logarithm of w(T) and removing the constant we get 



In I nq 



J-(S+l)/2 



1 



,j-(S+l)/2 



Kq- 



In Kq 



J-(S+3)/2 



j-(S+3)/2 



Kq 



This is a Ricmannian sum for an integral, and as e — > this yields, up to 
second order terms, 



1 



91n( K q- s / 2 +j 

h(x,t)— - 



hexagon 



where x, t are normalized coordinates inside the hexagon, h is the normalized 
height function, and j — x — t/2. 

Following [CKP] and [KOS] we know that the number of stepped surfaces 
in an e-neighborhood of a given height function is asymptotically 



6XP / a(Vh) 



where a is the surface tension that can be expressed through the Lobachevski 
function. (Note that the signs in [CKP] and [KO !] are different, we follow 
[CKP].) 

Consequently, the probability of the stepped surfaces in the e-neighborhood 
of a given height function is asymptotically proportional to 



exp 



ify a(vh)+ j h( X ,t) 



91n(/iq 



-S/2+j 



K q-S/2+i. 



and to find the limit shape one has to maximize this expression. We conclude 
that the limit shape can be found as a solution of the variational problem: 



o-(Vh) + / h(x,t) 



91n(Kq 



-S/2+j 



f<q 



-S/2+j 



di 
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Let us write down the Euler-Lagrange equation for this variational problem. 
It is convenient to use our parameter z instead of partial derivatives of the limit 
height function h. We know (see [KO]) that the Euler-Lagrange equation for 
the first term is the complex Burgers equation 

z 1 — z 

Adding it to the trivial Euler-Lagrange equation for the second term we 
obtain 

z 1 — z d] 

or 

Zt Z X K 2 q -S+2x- t + 1 

Z l-Z 11 K 2q-S+2K-t _ i 

This is a quasilinear equation. The equations for characteristics have the 
form 

dt _ 1 dx _ 1 dz _ K 2 q- s+2a; -' + 1 

ds z' ds Z— 1' ds ^ K 2 q-S+2a;-t _ i 

We find two first integrals 

zq l - K 2 q- s+2x (1 - z)q x 

u == V == . 

1 — zK 2 q~ s+2x-t ' 1 — ZK 2 q~ s + 2x ~ l 

Any solution of the partial derivative equation has the form (cf. the proof 
of Corollary 1 in [KO]) 

Q(u,v)=0, 

where Q is a suitable analytic function. This equation defines z as a function 
of t and x. 

In the g-Hahn case, it is known that Q(u, v) is a second degree polynomial. 
It is natural to assume that the same is true for the g-Racah case, and then the 
requirement that the (t, x)-curve where z degenerates to M is tangent to the six 
sides of the hexagon leads to the formula of Theorem 2.1. It remains a challenge 
to check if Q is still algebraic for more general polygonal domains, as was shown 
in [KO] for the g-Hahn case. 



3 Weight sums 

A horizontal lozenge in the tiling interpretation corresponds to a hole (^absence 
of a particle) in the nonintersecting paths or Appoint configuration interpreta- 
tion. The coordinates of the horizontal lozenge in the formulas (1) and (3) 
correspond to the coordinates (t, x) of a hole in the following way: 

i — t, 

j = x - t/2+\. 
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Now consider any vertical section t = to of the family of nonintersecting 
paths (see Figure 3) and fix the corresponding iV-point configuration X\ < x 2 < 
■ ■ ■ < x N . 

Denote by C(t; X\, . . . , x^) the product of weights corresponding to the holes 
on this vertical line. 

Denote by L(t; x\, . . . , xn) the sum of the products of weights corresponding 
to holes situated to the left of the vertical line. The sum is taken over all 
families of paths connecting the points {(0, 0), (0, 1), . . . , (0, N— 1)} to the points 
{(to,xi), (t ,x 2 ), (t ,x N )}. 

Denote by R(t; x\, . . . ,Xn) the sum of the products of weights corresponding 
to holes situated to the right of the vertical line. The sum is taken over all 
families of paths connecting the points {(to, (to, X2), ■ ■ ■ , (to, £j\r)} to the 
points {(T, S), (T, S + 1), . . . , (T, S + N - 1)}. 

The following three propositions correspond to the case of the g-Racah 
weight (3). 

Set 

K,s(x) = q- x + K 2 q x ~ S ~ t+1 
Here and below we use the g-Pochhammer symbol: 

(a;q) n = (1 - a)(l - aq) . . . (1 ~ aq^ 1 ). 

Proposition 3.1. We have 

L t (x\,...,x N ) =const t - } ] (fj, t ,s(xi) - Vt,s(xj)) 

\<i<j<N 

X 2 (q-^q^) t+N ^ Xi ■ (q;q) Xi ■ ^-t-S+l. q)t+N " 

Proof. Here and in the proof of the next lemma we should consider four cases 
depending on the value of t (see formulas (8)-(ll)). The proofs are similar in 
all four cases and we consider only the one that corresponds to the pictures in 
Lemma 10.2 and Lemma 10.3 (i.e., S < t < T - S) 

Let us use Lemma 10.2 that expresses the required weight sum in terms of 
the point configuration complementary to {xi}, i.e., in terms of the positions 
of the horizontal lozenges that we call holes. Denote the positions of holes by 
{yi}. Performing the limit transition (2) we get. 



L t (xi,...,x N ) = const- Yl i^t,s{yi) ~ Vt,s(Vj)) 

l<i<j<S 

tt q {1 ,q)s+N-i- yi {q /« ) q)s+N-i- yi 

t <J< s (q;q)s+N-i-yA<i~ 2N+1 / K2 ;<i)s+N-i-y t 

To finish the proof we rewrite the last formula in terms of particles {xi} 
instead of holes {yi}. For the second factor this procedure is simple, while for 
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the first one we use the following observation: 

n (vt,s(yi) - tH,s(vj)) = n i^t,s(xi) - nt,s(xj)) 

1<i<i<N l<i<j<N 

x n (^t,s{ u ) - tH,s{ v )) 

0<u<v<S+N~l 

1 



n 



- - v 0<u<xi xt<u<S+N-l 

The product over u < v depends only on t, while the last two products over u 
are easily expressible in terms of g-Pochhammer symbols. 

□ 

Proposition 3.2. We have 

R t (xi,...,x N ) = cemstf JJ (fH,s( x i) ~ ^sfaj)) 

l<i<j<N 

N (1 - K 2 q 2 *i-S-t+i\qXi(T-t+N-i) 

X H {q^iq^s+N-i-xi ■ (q;q)x z +T-t-s ■ {K 2 q x -- T + 1 ;q) N+T -t ' 

Proof. Performing the limit transition (2) in Lemma 10.3 we get 
R t (xi, ... ,x N ) = const ■ ]~J {fH,s(Vi) ~ Mt,s(%)) 

l<i<j<S 



n 



g (S-T +t )(S+iV-l- W )(g-iV-S + l ;5)g+JV _ 1 _ s/i(g T-S-JV + l /K2 . 9)s+iV _ 1 _ i/4 



Again expressing all the factors in terms of {xi} we obtain the desired result. 

□ 

Proposition 3.3. We have 

N 



i=l 

Proof. Clearly, 

S 

-S/2-t/2+l/2 



C t (x u . . . , x N ) = constt • [] 1 _ K 2 g L-s-t+i • 



c t (xx, . . .,x N ) = Y\_ ( K q v% 



i=l 



KqVi -S/2-t/2+l/2 j ' 



where {y.i} is the point configuration complementary to {xi}. Expressing the 
product in terms of {xi} we get the result. □ 
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4 Distributions and transition probabilities 



We consider our probability measure on the set of tilings of a given hexagon as 
a Markov chain in the space of iV-tuples of integers. The Markov property can 
be easily seen in the following form: The past and the future are independent 
given the present. In this way the Markov property reduces to the fact that a 
lozenge tiling of the hexagon is a union of a tiling to the left of a vertical line 
t = const and a tiling to the right of this vertical line. Denote this Markov 
chain by X(t), t = 0, 1, . . . T. 
Set 

Xf}* T = {xeZ: max(0, t + S - T) < x < min(< + TV - 1, S + JV - 1)} 

and 

X N% = {{ X ll X 2, • ■ ■ ,X N ) € (X^ T ) N : X! < x 2 < ■ ■ ■ < x N }. 

St- 

Xjj ji is the section of our hexagon by the vertical line with coordinate t, and 
s t 

X N ' T is the set of all TV-tuples in this section. 
Clearly X(t) takes values in X^'tp. 

The following theorem gives the one-dimensional distributions of our chain at 
the time t, which are probability distributions on TV-tuples of integers. Denote 
by ps,t the one-dimensional distribution of the process X(t) (below we explain 
why we keep S in the notation but omit all the other parameters). 

Theorem 4.1. 

N 

Pmb{X(t) = (x 1 ,x 2 , ■ ■ -,x N )} = const ■ TJ(^ tiS (x 4 ) - (i t , s (xj)) 2 w t ,s(xi), 

i<j 2=1 

where 

IH,s(*) = <T* + *V~ 5 - t+1 

and 

(—iy+ s q x ( 2N + T - 1 )h — K 2 q 2x ~ t ~ s+1 ) 

Wt ' S ^ (q; q) x (q; q) T -s~t+x(q~ 1 ; q'^t+N-x-iiq- 1 ; q^ 1 ) s+N-x-i 

1 



(«V" T+1 ; ?)T+w-t(KV-*- s+1 ; q)N+ t ' 

Proof. Clearly, 

Prob{X(t) = (xi,x 2 ,-.-,x N )} oc L(t;xi, . . . , x N )C(t;xi, . . . ,x N )R(t;xi, . . . ,x N ) 

Propositions 3.1, 3.2, 3.3 imply the result. □ 

Observe that w t .s(x) is (up to the factor not depending on x) the weight 
function of the g-Racah orthogonal polynomials, see e. g. [KS, Section 3.2]. 
These polynomials are given by the formula 

RnMxba, 0,1,8 | q) = 4^3 ( q n,a %^ X iq lSqX+1 \ , (7) 



14 



where 



fj,(x) 



-fdq 



x+l 



and aq = q M or j38q = q M or jq = q M for a nonnegative integer M. They 
are orthogonal on {0, 1, ... , M} with respect to the weight function 



w(x) 



(aq,P6q,^q,-fSq; q) 3 



1 - 1 5q 2x+1 



(q, a 1 j5q, fi-^q, Sq; q) x (a(3q) x (l - jdq) ' 



The correspondence between the parameters of polynomials and parameters 
of our model is established in the following way. 

(i). t<S, T - t - S > 0, 0<x<t + N-l, 



Q(qRacah) 
(X(qRacah) 
P(qRacah) 
l(qRacah) 
^(qRacah) 



-S-N 
^S-T-N 



K 2 q 



-N 

-S+N 



(8) 



(ii). S-l<t<T-S + l, 0<x<S + N-l, 



Q(qRacah) 
^■(qRacah) 
P(qRacah) 
'y(qRacah) 
^(qRacah) 



q -t-N 
q t-T-N 
q -S-N 
K 2 q -t+N 



(9) 



(iii). T - S - 1< t < S, 0<x-(t+S-T)<T-S + N-l, 



Q(qRacah) 
Ot{ q Racah) 
P(qRacah) 
1(qRacah) 
S(qRacah) 
% (qRacah) 



= q 



-T+t-N 
-t-N 



T-N+S 



K 2 q -T+t+N 



(10) 



t- S + x 



(iv). t > T - S - 1, t>S-l, 0<x-{t + S-T)<T-t + N-l, 



Q(qRacah) 
^(qRacah) 
P(qRacah) 
1(qRacah) 
d(qRacah) 
%(qRacah) 



= q 



-T-N+S 

-S-N 

-T+t-N 



(11) 



= K 2 q~ T + N + s 
= T-t- S+x 
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Let us also describe what happens if one performs the limit transitions de- 
scribed in Section 2.2. 

If we send k —> (the weight of a plane partition becomes proportional to 
q-voiume^ t nen we obtain the weight function 

qX(2N+T-l) 

(q; q)x{q; q)T-s-t+x(q- 1 ;q- 1 )t+N-x-i(q~ 1 ;q~ 1 )s+N-x-i 

This is exactly (up to the factor not depending on x) the weight function of the 
g-Hahn polynomials. 

If we set k = q K and send q —> 1, then the weight function becomes 

w{x) 



x\(T -S-t + x)\{t + N-x- 1)!(5 + N - x - 1)! 

K + 2x - t - S + 1 



(K + X-T+ l) T+N -t(K + x-t-S+ l) N+t ' 

This is the weight function of the Racah orthogonal polynomials. 

Finally if we send n — > and set q — 1 (this case corresponds to the uniform 
measure on the set of lozenge tilings of the hexagon) then we get the weight 
function 

= x\(T-S-t + x)\(t + N-x-l)\(S + N-x-l)\' 

This is the weight function of the Halm polynomials. This case was previously 
studied in [Jl], [J2], [Gor], see also references therein. 

We also need the transition probabilities of the Markov chain X(t). 
Proposition 4.2. 

PTob{X(t + 1) = Y\X(t) = X} 



where 



, TT ft+l,SWJ — h l t+l,S\!Jj ) TT / \ TT / n 

co ^-n , tAxi) -, tAX]) n n «*m> 

1 _ .2 s+JV-f 

Wo (x) = -(l-^ +T - t - S ): 



1 - ft^s-i-SH-l 

and 



_ -2„x-T+l 

Wl (a;) = g T+jV - 1 - t (l- (Z a; - s - JV+1 ) ' ' 



1 _ K 2^2K-t-S+l ' 



Proof. We use 

Lt(X)C t (X)C t+1 (F)ii t+1 (y) 



Prob{X(t + 1) = F|X(i) = X} 



L t (X)C t (X)J? t (X) 

R t+1 (Y)C t+1 (Y) 



Rt{X) 

and Propositions 3.3, 3.2. □ 
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Next let us compute the cotransition probabilities (t — > t — 1). 
Proposition 4.3. 

Vvob{X{t - 1) = Y\X(t) = X} = 

T-r fH-i,s(Vi) - IH-\,s{yj) ~ , \ tt ~ / \ 
const- — — - — Wi{Xi) w (Xi), 



where 

tD (a:) = -(l-« s - t - iV+1 ) 



1 _ K 2 n x-S-t+l 
x-t-N+l\ 1 K 1 



1 - K 2 q 2x-t-S+l 

and 

m(x) = q- { - t+N - x \\-(f)-. 



1-K 2 q x+N - S 



1 - K 2 q 2x-t-S+l 

Proof. We use 

Fl ° h{X(t - 1} = =X}= L t (X)C t (X)R t (X) 

Lt{X) 

and Propositions 3.1, 3.3. □ 

5 Families of stochastic matrices 

This section and the next one are similar to [BG], where the Hahn case was 
treated, and we have tried to keep the notations and statements of theorems 
unchanged where possible. 

5.1 Definition of matrices 

St St St St 

We want to introduce four families of stochastic matrices P t + , P t l , P s ' + , Pgl- 
Pf_f(X,Y) is an x \X S ' t+1 \ matrix, X = {x x < ■ ■ ■ < x N ) G X s,t , 

Y = ( Vl < ■ ■ ■ < y N ) G X s ' t+X ; 

if yi — Xi G {0, 1} for every i, then 

P ti (X,Y)^ const- [[————— [[ Wl { Xl ) [[ w ( Xi ), 

iKj «.fF»J Vt^^]) y i=Xi+1 y. =X( 



where 



i _ K 2 n x+N-t 

w a (x) = -(l- q * +T - t - s )-_ 1 ~- ' 



I _ K 2g2x-t-S+l ' 

1 _ K 2 x-T+1 
Wi(x) = q^- 1 -^! - ga-S-AH- 1 ): 1 " 1 



I _ K 2g2x-t-S+l ' 
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and P t S +(X, Y) = otherwise. 

Pf;*(X,Y) is an x \X S+1 ^\ matrix, X = ( Xl < ■ ■ ■ < x N ) G A ,S '*, 

r=(y 1 <---<y„)GA' s + 1 ^; 

If Vi — Xi £ {0, 1} for every i, then 

Mt,s+i(y») - in,s+i{yj) 



p§>*{x,y) = const -n ^ s+i * - n n 



where 

1 - K 2 ^ +Ar - S 



u; (a:) = -(1 - <f+ T -*- s ) 



I _ K 2g2x-t-S+l ' 



and Pg'*(X, Y) = otherwise. 

P t ^*(X,y) is an lA" 5 '*! x I* 5 '*- 1 ! matrix, X = (an < • • • < jcjv) G A 5 ' 1 

r = (y 1 <---<y„)eA'^- 1 ; 

If yi — Xi £ {—1,0} for every i, then 

pff ( x,y) = const . n " Mt T ( f } n fiiNnw. 



where 



i<j ' v ' \ J / y i=Xi — l Vi=Xi 



i _ J2 n x— S— 1+1 

*o(x) = -(l-^-*- Ar+1 ) i " " 7 



I _ K 2g2x-t-S+l ' 

1 - K 2 q x+N - S 



and P^pf, y) = otherwise. 

p£_*(X,y) is an \X S >*\ x lAf 5 " 1 -*! matrix, X = ( Xl < ■ ■ ■ < x N ) G # s >' 
K = (yi < ■ ■ ■ < y„) G A" 5 ^ 1,4 ; 

If iji — Xi G { — 1, 0} for every i, then 

pi±{x, y) = const . n ^-^i " ^-y y n fiiNlI^ 



where 



1 _ K 2 z-S-t+1 



U>1 



(.x) = g -( s+Ar - 1 )(l- g a; ) 



I _ K 2g2x-t-S+l 
I _ K 2qX+N-t 



I _ K 2g2x-t-S+l ' 

and Pg'*(X, Y) = otherwise. 

Looking at the sets that parameterize rows and columns of these matrices 

St St St 

one can say that P t + increases t, P t l decreases t, while P s ^ increases S and 

S t 

P s !_ decreases S. This explains our notation. 
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Theorem 5.1. With appropriate choices of normalizing constants, all four types 
of matrices defined above are stochastic. They preserve the family of measures 
ps,t- In other words 

P?± t (X,Y) = l, J2 Pf± t (X,Y) = h (12) 



Ps,t±i(Y)= Pt±\X,Y)- ps,t{X), 

Ps±iAY)= £ Psi(X,Y)-p s AX). 

Proof. Propositions 4.2 and 4.3 imply the claim for P t +(X, Y) and P t J t (X, Y). 

Now observe that the space X s ' 1 is unaffected when we interchange param- 
eters t and S, i.e., 

Moreover, the measures ps t t are also invariant under S <-> i, i.e., 

Ps.t = Pt,s- 

(This is a consequence of our special choice of the parameter k which included 
additional factor q~ s / 2 .) 

Finally, note that Pf+(X, Y) becomes Pff (X, Y) under S <-> t and P t ?l* (X, Y) 

becomes Pfl*(X, F). 

Therefore, applying 5 <-> i to the relations for P t ± we obtain the needed 
relations for Ps±. □ 

5.2 Determinantal representation 

In this section we write our stochastic matrices in a determinantal form. This 
representation is very convenient for various computations. 

First, we introduce 4 new two-diagonal matrices. 

For x € X s -*, y G X S ' t+1 , 

q T + N ^-\l qS-N+l) , ify = x+l, 

U&\x, y) = \{l- ( f +T - t - s ) 1 ^X£^ , i£y = x, 

otherwise; 




0, otherwise 



T-t-S\ 1-K 2 q* " 
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for x G X 5 '*, y G X s '*" 1 , 

2 x + JV-S 



U t S l\x,y)={ 



k 0, otherwise; 
and for a: G X 5 '*, y G X s - M , 

r g -(5+iv-i) (1 _ g x ) _l^^ r ; if y = a:-l, 

Ut(x, y) = l-(l- f-s-N+i) }:$fc_ s -Z , ify = x , 

[0, otherwise. 

St St 

It is possible to express the stochastic matrices P t ± , Pg ± as certain minors 
of the matrices defined above. 

Proposition 5.2. We have 

P t + (X, y) = const ■ M ^— ^ ^-det[£/ t+ (aii.i/j-) ij=i,...,iv 

vt,s{xi) - m,s{xj) 

nS.tr v . TT Mt,S+l(2/i) — Mt,S+l(?/j) , ,r TT S.tr m 

«<j ^t^sy^i) H-tjSy-Lj) 

u S,tr v v \ . TT ^t-l,s{Vi) — fM-l,s{yj) A ,r TT S,tr 

*t- (*> r ) = cons< ' 11 r^T~ det[[/ t i I/j) ij=i,...,iv 

r>S.tr v ^r\ . TT Ht,S-\{yi) — Mt,S-l(%) , ,r T TS,tr \i 

P s l(X,y) = const- — ^ / dct[L/g_(xi,yj)]»,j^i....,Ar 

f<* tH,s{Xi) - IH,s(Xj) 

Proof. Straightforward computation using the definitions of the stochastic ma- 
trices Pf±, Pg± and the matrices Uf±, Ug±. 

Any submatrix of a two-diagonal matrix, which has a nonzero determinant, 
is block-diagonal, where each block is either an upper or a lower triangular 
matrix. Thus, any nonzero minor is a product of suitable matrix elements. □ 



5.3 Commutativity 

St St 

Theorem 5.3. The families of stochastic matrices P t ± and P s ' ± commute, that 
is 



p s,t 


pS,t+l 

r s- 


_ p S,t 


pS-l,t 


P s,t 


pS,t-l 


_ p S,t 


pS-l,t 


pS.t 


pS.t+1 


_ p S,t 


pS+l,t 

r t+ 


pS,t 

r t- 


pS,t-l 

^s+ 


pS,t 

- ^s+ 


pS+M 
r t- 



for any meaningful values of S and t. 
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Proof. Proofs of all four cases are very similar and we consider only the first 
one. 

{p s,t . (X) Y)= £ p st z) . p s,t+i (Z) y) 

zeArs,t+i 

x det[t/ t S | :*(x l ,z J )] 4J= i^.. iA rdet[?7|^ +:L (z 4 ,y i )] ii: , = i i ... jA r. 

Applying the Cauchy-Binet identity we obtain 

det[U t + {xi, Zj)]ij =1 ^.. !N det[U s '_ (^,j/j)]i,j=i,...,jv 

= det[(^ 4 • C/|l t+1 )(x 4 , % )] JJ=1 ,..., w . 

Thus, 

. TT Mt+l,S-l(j/») _ ^t+l,S-l{Vj) , , rm S,t rr S,t+lw \i 

= const- 1| /v^_„ ^ det[(tf t+ ■ C/^l )(a; i ,y j -)]ij=i,...,jv. 

Similarly, 

(p |,. p 5 + -M )(X;r) 

. TT IM+l,S-l(yi) - IM+l,S-l(yj) , , UTT S,t TT S-l,u, \n 
= const || — — — — det[(U s _ ■ U t+ ){x h yj)]ij=i,...,N- 

Our claim reduces to verifying the equality 

£# . u§± +1 = uj£ ■ d*- 1 -*. 

Note that this will also imply the coincidence of normalization constants, since 
all matrices under consideration are stochastic. 
A straightforward computation yields 

rrS,t j-j-S,t Tj-S.t-\-l jjS.t TjS — 1}^ 

u t+s- ~ u t+ ■ u s- — u s- ' u t+ ■> 

where 

' ui, ify = £ + l, 



u , ify = x, 
u-i, ify = x-l, 
0, otherwise, 
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and 



"1 



JT+N-l-t 



(1-q 



x-S-N+l 



)(i-g 



x-S-N+2\ 



(1 



n 2 q x ~ 



S-t+i 



)(1 



(1 



u = -(l-q x 



-S-N+l 



i - kV+ jv -« 

■ 2g2x 
1 - K 2 ^- T+1 



1 - K 2 g 2x-t-S+l 
■'- s "'(I - ,/ - 1 ) J ~ + C 1 _ g a+T ~*~ 5 ) 1 1 ~ <g 2 to-/ 



□ 



6 Perfect sampling algorithm 

6.1 Definition of transition matrices 

In this section we aim to define two new stochastic matrices 

PLs+i(X,Y), XeQ(N,T,S), Yen(N,T,S+l) 

and 

P§„s-i{X,Y), XeSl(N,T,S), YeSl(N,T,S-l) 

that preserve the measures n(N,T, S,q, n). Both P§^ s+1 and P§^g_ 1 depend 
on parameters N, T, q, k but we omit these parameters from the notation. 

Suppose we are given a sequence X = (X(0), X(l), . . . , X(T)) £ Sl(N,T, 5) 
(recall that X(t) £ X s,t ). Below we construct a random sequence Y = (5^(0), . . . , Y(T)) £ 
Sl(N,T, 5 + 1) and therefore define the transition probability (or, equivalcntly, 
stochastic matrix) Pg l _ JrS+1 (X,Y). 

First note that F(0) £ X s+1 '° and \X S+1 >°\ = 1. Thus, Y(0) is uniquely 
defined. We will perform a sequential update. Suppose Y(0), Y{1), . . . , Y(t) 
have been already defined. Define the conditional distribution of Y(t+ 1) given 
X, Y(0),Y(l),...,Y(t) by 



Prob{F(i + 1) = Z} 



P f+^(Y(t),Z) ■ P§+ 1 ' t+1 (Z,X(t + l)) 
(p S ++Mp |+M+i )(r(<);X(t+1)) 



P|; t+1 (X(i + 1), Z) ■ F t 5+1 ' t+1 (Z, Y(t)) 
(P^ +1 P t S + 1A+1 )(X(t+l),Y(t)) 



~ ' (13) 

(The second equality follows from p s+1 ,t + i(X)Pf+ 1,t+1 (X, Y) = p s+ltt (Y)P t s + ht (Y, X).) 



This definition follows the idea of [DF, Section 2.3], see also [BF]. 



22 



Observe that (p£ +M pf + M+1 )(Y(i), X(t + 1)) > (here and below see [BG] 
for more details). 

One could say that we choose Y(t + 1) using conditional distribution of the 



middle point in the successive application of P t ^_ +1 '* and p^ 1 ' 1 ^ 1 ( or p s+ 
and p t ^ 1 ' t + 1 ) ; provided that we start at Y(t) and finish at X(t + 1) (or start 
at X{t + 1) and finish at Y(t)). 

After performing T updates we obtain the sequence Y. 
Equivalently, define P§^ s+1 by 

fr-i P t s + ht (Y(t),Y(t + 1)) • Pf + ll4+1 (Y(f + l),X(t + 1)) 



,S+l,t+l 



P§„ S+1 (X,Y) = { 



(Pi 



S+l,t pS+l,t+l 



T-l 



S- 



)(F(t),X(t+l)) 



if II {P^ 1 ' t P S s + X ' t+X ){Y{t),X{t+ 1)) > 0, 



t=o 

0, otherwise. 



Theorem 6.1. T/ie mate Pf^ s+1 on fi(JV, T, S) x Q(jV, T, 5+1) is stochastic. 
The transition probabilities P§,_^g +1 (X, Y) preserve the measures n(N, T, S, q, k): 

»(N,T,S+l,q,K)(Y)= ]T P§„s+i(X>Y)n(N,T,S,q,K)(X). 

xen(N,T,s) 



Proof. See [BG]. 

Similarly to Pg^s+i, one defines a transition matrix 



□ 



P^ S -i(X,Y), XeSl(N,T,S), Y etl(N,T,S-l), 



by 



(T-i if-^fy^Yit + 1)) ■ Pl- x ' t+1 (Y(t + \).X(t + 1)) 



n 



S-l,t pS-M+1 
t+ +S+ 



)(Y(t),X(t + l)) 



T-l 



if n {p t s + M Pi+ M+1 )(n*).^+i)) >q, 

t=0 



0, otherwise. 

Similarly to (13) there is another way to write P§^s-i because of the equality 



P 



^(myfi + i)).^ 



s-i,t+i(Y( t + l),X(t + l)) 



(Pi 



S-l,t p S-l,t+l 

J- I 



)(Y(i),X(i + l)) 



P, 



S,t+1 



s- 



(x(t + i),Y(t + i))-p; 



S-l,t+l 



(Y(t+l),Y(t)) 



^pS,t+i p s-i,t+i ){x{t + 1 ^ Y{t)) 



Similarly to Theorem 6.1 one proves the following claim. 
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Theorem 6.2. The matrix P§„ S -i on n ( N > T > S ) x fi W T > ^ *) * s stochastic. 
The transition probabilities Pg^g_ 1 (X, Y) preserve the measures n(N, T, S, q, n): 

n{N,T,S-l,q,K)(Y)= PLs-i(X,Y)^N,T,S,q }K )(X). 

xen{N,T,s) 

Remark. The above construction performs a sequential update from t = to 
t = T . One can equally well update from t = T to t = by suitably modi- 
fying the definitions. The resulting Markov chains also preserve the measures 
H(N, T, S, q, k), and they are different from the Markov chains defined above. 

6.2 Algorithm for the S i-> S + 1 step. 

Now we suggest an algorithmic description of the Markov chain from the previ- 
ous section. 
Denote 

J T f „ K C T"| = 1 g l-k q 1 g 

^r-t-s-i^ _ i _ K 2 qX -T+i 1 _ K 2 q 2 X -ts-i 

and 

k 

P(x, t, q, k, S, T; k) = \\_p{x + i — l,t,q, K, S,T) 

i=l 

{qX+ T-t-S-l. q)k {K 2 qX -S-t-l. q)k ( K 2 q 2 X -t-S+l. q 2 )k 
q k(T-t-S-l)( q x+l. q)k ( K 2 qX -T+l. q)k ( K 2 q 2 x -t-S-l. q 2 )k - 

Denote by D(x, t, S; n) (it also depends on q, k, T, but we omit these param- 
eters) the probability distribution on {0, 1, . . . , n} given by 

(14) 

Suppose we are given X = (X(0),X(1), X{T)) e Q(N, T, S). We want 
to construct Y = (T(0), Y(l), Y{T)) e Q(N, T,S+ 1). 
In the first place we note that Y(0) is uniquely defined, 

F(0) = (0,1,...,JV-1). 

Then we perform T sequential updates, i.e., for t = 0, 1, . . . T — 1 we construct 
Y(t + 1) using Y(t) and X(t + 1). Let us describe each step. 

Let Y(t) = (j/i < y 2 < ■ ■ ■ < Un) and X(t + 1) = (xi < x 2 < ■ ■ ■ < x N ). We 
are going to construct Y(t + 1) = (z\ < z 2 < ■ ■ ■ < z N ). 

Recall that 

zt e X 5+M+1 = {x e Z | max(0, t + S-T + 2)<x< min(i + N, S + N)}. 
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Observe that Y(t) and X(t+1) satisfy (p£ +M Pf + M+1 )(T(>), X(t+1)) > 0. 
This implies that Xi — yi is equal to cither —1, or 1 for every i. 

• First, consider all indices i such that x% — yi = 1. For every such i we set 

• Second, consider all indices i such that Xi — yi — — 1 and set Zi = yi- 

• Finally, consider all remaining indices, i.e., all i such that Xi = yi. Divide 
the corresponding Xi's into blocks of neighboring integers of distance at least 
one from each other. Call such a block a (fc, Z)-block, where k is the smallest 
number in the block and I is its size. Thus, we have 



Xi = yi = k, x i+ i = y t+ i = k + 1, 



Xi+i-i = y l+ i-i =k+l-l 



and 



l < k - 1, > k + l. 



For each (£;, Z)-block we perform the following procedure: consider a random 
variable £ distributed according to Z?(fc, 5; I) (£'s corresponding to different 
(fc,/)-blocks are independent). Set = Xi for the first £ integers of the block 
(their coordinates are k,k + 1, . . . ,k + £ — 1) and set z, = + 1 for the rest of 
the block. 

At Figure 4 we provide an example of constructing Y(t + 1) using X(t + 1) 
and Y(t): there is only one (k, /)-block and it splits into two groups, here £ = 2. 



X(t) X(t+1) 



Y(t) Y(t+1) 



o o 



o o 



o 

i 




could not jump (first case) 



split point determined by (8) 



(block case) 



~ was forced to jump (second case) 



o o 

S+l 



Figure 4. Example of (fc, /)-block split, I = 4, £ = 2. 



Theorem 6.3. The algorithm described above is precisely the S i— > S+l Markov 
step given by PjLvS+i- 



Proof. Straightforward computations. See [BG] for some details. 



□ 
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Remarks. Setting n = in the formulas for the distribution D(x,t, S;n) 
we obtain the perfect sampling algorithm for boxed plane partitions distributed 

no volume 
as q 

Sending q — > 1 in the formulas for the distribution D(x, t, S; n) as described 
in Section 2.2, we get a perfect sampling algorithm for the Racah case (recall 
that in this case the weight of a horizontal lozenge is proportional to a linear 
function of its vertical coordinate). 

6.3 Algorithm for S h- > S — 1 step 

Using similar methods we can also obtain S S — I Markov step which gives 
alternative way to sample a random tiling: We start from the case T = S and 
then perform some amount of S i— > S — 1 steps. 

The step algorithm is very similar to the S > S + 1 one. 

Denote 

qt+i-Sfj _ qX-t-N-i\ j _ K 2qX+N-t-i ^ _ K 2 q 2x-t-s+i 

p(x, t, q, K, S, T, N) — (1 — gX-S-N+lj i _ K 2 q x+N-S+1 \ _ K 2 q 2x-t~S~l 

and 

k 



P(x, t, q, k, 5, T, N; k) = ~^^p{x + i — l,t,q,K,S,T, N) 



_ g fc(*+l-S)(^-t-jV-l ;g ) fc ( K 2 ga +iV-t-l ;g ) fc ( K 2 g 2tt-t-S+l ;g 2) fc 
(qxS-N+l. q)k {K 2 q x+N-S+l. q)k {K 2 q 2x-t-S-l ;q 2 )k ■ 

Denote by D(x, t, S; n) the probability distribution on {0, 1, . . . , n} given by 

Prob({fc}) = D(x, t, S; n){k} = P^^^S^N^) ^ 

Ej=o p ( x i *5 1, «j 5", T, AT; j) 

Suppose we are given X = (X(0), . . . , X{T)) G n(AT, T, 5). We want 

to construct Y = (F(0), . . . , F(T)) G fl(N, T,S— 1). 
As above, note that F(0) is uniquely defined, 

F(Q) = (0,l,...,iV-l). 

Then we again perform T sequential updates, i.e., for t = 0,1,... T — 1 we 
construct Y(t + 1) using Y(t) and X(t + 1). Let us describe each step. 

Let Y(t) = (yi <V2 <■■■ < Vn) and X(t + 1) = (xi < x 2 < ■ ■ ■ < x N ). We 
are going to construct Y(t + 1) = [z\ < z 2 < • • • < zjy). 

Recall that 

Zi G = {x G Z | max(0, i + S - T) < x < min(f + AT, S + N - 2)}. 

Y(t) and satisfy (P t + _M Pf+ M+1 )(F(*), X(t+1)) > 0. This implies 

that Xi — yi is equal to either 0, 1 or 2 for every i. 
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• First, consider all indices i such that Xi — jji = 0. For every such i we set 

Zi Xi . 

• Second, consider all indices i such that x\ — = 2 and set z% = y%-\-\. 

• Finally, consider all remaining indices, i.e., all i such that Xj, = yi + 1. 
Divide the corresponding Xj's into blocks of neighboring integers of distance 
at least one from each other. Call such a block a (k, Z)'-block, where fc is the 
smallest number in the block and I is its size. Thus, we have 

Xi =yi + l=k, x i+ i = y l+ i + 1 = k + 1, . . . , Xi+i-i = Vi+i-i =k + l-l. 

For each (fc, Z)'-block we perform the following procedure: consider random 
variable £ distributed according to D(k, t, S; I) (£'s corresponding to different 
(k, I) '-blocks are independent). Set Zi — yi for the first £ integers of the block 
(their coordinates are k — l,fc,...,fc + £ — 2) and set z< = yi + 1 for the rest of 
the block. 

Theorem 6.4. The algorithm described above is precisely S i— > 5 — 1 Markov 
step defined by Pg 1 _ > g_ 1 . 

The proof is similar to Theorem 6.3. 
6.4 Markov evolution of the top path 

The S 1 1—* S + 1 Markov step described in the previous section has the following 
property: Its projection to the set of topmost horizontal lozenges (or the top- 
most holes in terms of nonintersecting paths and point configurations) is also a 
Markov chain. This Markov chain is an exclusion type process. Let us describe 
it. 

The general setting is as follows. The state space of our discrete time Markov 
chain consists of semi-infinite particle configurations {ei < ei < < . . . } in 
Z. At each time moment every particle either stays or jumps to the left (any 
distance) avoiding collisions and jumps over neighbors. Jumps are performed 
sequentially. First, the leftmost particle (ei) jumps, then the second one and 
so on. The distribution D of the length of the jump of a particle depends on 
the number of the particle, moment of time, current position of the particle (e,) 
and the distance between the current position of the particle and the position 
of the previous particle (ej_i) in the next moment of time. At time we have 
the step initial condition, i.e., ej = i + const. 

Now let us turn back to our situation. All particles are enumerated by the 
parameter t and our time parameter is S that changes from to T. Consider 
a sequence {uf }t=i .., where uf is the vertical coordinate (in our notation - 
x) corresponding to the topmost hole inside the hexagon for t < S and uf = 
N+t—1 for t > S. (We can also view uf as the vertical coordinate corresponding 
to the tth hole, if we count all holes, not just the ones inside a hexagon, starting 
from the line x — 0.) 

The evolution of {uf } is precisely our Markov process. When S — the 
configuration consists of points N,N + 1, N + 2, . . . . The distribution of the 
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length of jump of the particle with coordinate uf at the time moment S is given 
by the distribution D(wfj t 1 1 + 1, t, S; uf - uf+l - 1) (see (14) for the definition). 
Note that when S = T the configuration consists of points 0,1,2,.... 

We can also obtain in a similar way a Markov chain for the bottommost 
holes. 

Finally, we may construct two more similar processes using the S i— » S — 
1 Markov chain (for this chains the direction of particle jumps changes and 
distributions D are replaced by distributions D). 

7 Correlation kernel 

The aim of this section is to obtain the formulas for the correlation functions 
of random point configurations in 1? obtained from the random tilings we are 
interested in. 

7.1 Expression via orthogonal polynomials 

Recall that a tiling of a hexagon corresponds to some family of nonintersccting 
paths that can be viewed as a point configuration in Z 2 . Let us denote this 
configuration by M. 




o 

o o o 
o o o o o 
o » o ° o 
o o o o 
o o 



As above, we denote the horizontal coordinate by t and the vertical coordi- 
nate by x. 

We want to compute the correlation functions of this random point config- 
uration. 

Recall that the nth correlation function is defined by 

p n (t 1 ,x 1 ;...;t n ,x n ) =Prob{(ti,Xi) G M, . . . , (t„, x n ) G M} 

for any collection {(£j, a^)}t=i,. „ of distinct points in Z 2 . 

To compute the correlation functions p n we are going to use a variant of the 
Eynard-Mehta theorem (see [EM] and [BO, Section 7.4]). Let us state it first. 

Proposition 7.1. Assume that for every time moment t we are given an or- 
thonormal system {/* } n >o i n 'a({0, 1, • • • , L}) and a set of numbers Cq, c\ , . . . . 
Denote 

v t , t+1 {x,y) = Y J C t J t n {x)f t n +1 {y). 

n>0 
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Assume also that we are given a discrete time Markov process Vt taking 
values in N -tuples of elements of the set {0, 1, . . . , L}, with one-dimensional 
distributions 

(det [/tife)] iJ=1) ... >iV W 

and transition probabilities 



det [v t)t+l {x iiy j)] i j=l _ N det [ft-l{Vj)] itj=1 

JV-l 

det^.!^)] 11 ci 



n=0 

Then 

Prob{xi € Vki,- ■■■,x n £ V kn } = det [K(k l ,x l ; kj,Xj)] iJ=h 

where 

N-l 

K(k,x;l,y)=J2 -i*fi(x)fi(v),k>l; 



K(k,x;l,y) = -J2 <£'7i*(*)/i(v)> k < l > 

i>N 

^* = 1, ^' = cf ■ cT 1 c\-\ 

Theorem 7.2. The Markov process X(t) meets the assumptions of Proposition 
7.1. 

The Markov process Vt is precisely our Markov process X(t). The orthonor- 
mal functions f^(x) are the normalized g-Racah polynomials multiplied by the 
square root of their weight function (see Section 4 for the definition of g-Racah 
polynomials, their weight function, and the correspondence between parameters 
of these polynomials and our parameters t, q, n, N, T, S): 

where (R^, i£„) is the squared norm of the q-Racah polynomials with respect to 
the weight function wt,s(x). This norm can be obtained from the norm of the 
g-Racah polynomials provided in [KS] (w tt s( x ) differs from the weight function 
of [KS] by a factor not depending on x). The explicit formula is a little bit 
different in the four cases of correspondence between parameters of polynomials 
and t, q, k, N, T, S. For instance, in the case given by formula (8): 

(R t R t , (-ir s (l- 2N - T+2 ,k-* q S- N ; q) t+N ^ 

[ nl { K - 2 q- 2N + 1 ,q s - T - N + 1 ,K 2 q-t- s + 2 ,q- t - N +l; q) t +N-i 

(1 - q- T - 2N )(q, q- T - N + t +\ K - 2 q- 2N + 1 ,q s - T - N ; q) n 
X (1 - q-2N-T+2n+l^ q -S-N^ q -2N-T+l^ K 2 g -T q - t -N+l. q ) n 

K 2n q -n(S+t+l) 



(q; q)T-s-t(q~ s ~ N+1 ;q)s+N-i(K 2 q - T+1 ;q) T +N-t ' 
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However, this long formula is not important for us, since factors involving 
it always cancel out. In particular, in the case given by (8) the quotient 
R t ^ rl )/(R t n , Rn), that is crucial for us, is simply 

(R t +\R t + 1 ) _ (I - qT+N-t-n-l^ _ q -t-N +n) 

(Ri,Ri) ~ {l-q-*-»)* 
The constants c* are given by 



4, = - - qT+N-t-n-iy {U) 

Proof. Theorem 4.1 yields 

1 n 

On the other hand 

jv 

det [fi-i{xj)] itj = h ... >N = const ■ Y[\Jw t ,s( x i) det \ R ^ x i\ =i,...,N ■ 

i=l 

The last determinant is a Vandermonde determinant in variables a t g{xj), hence 

2 1 N 

(det [fi-i( x j)] itj =i t .^ N ) = zjJJ(l*t,s( x i) ~ Vt,s{ x i)Yn w t,s{xi)- 

i<j i—1 

Coincidence of the constants (Z — Z') follows from the fact that the left-hand 
side in the last equality defines a probability distribution. 

Thus, the one-dimensional distributions of our process have the required 
form. 

Next, we need the following standard facts. 
Lemma 7.3. The following relation for the basic hypergeometric function holds: 

(c- w)(l - d)i<f>3 



a, b, c, qd 
u, v, qw 



q;q) +{w- d)(l - c) 4 03 
(c- d)(l - w) 4 (/) 3 



a, b, qc, d 
u, v, qw 
a, b, c, d 
u, v, w 

In terms of q-Racah polynomials, this relation can be rewritten as 



q;q 

q;q). (18) 



(q x - 97) (1 - 7Sq x+1 )R n U(x); a, /3, qj,8\q 

+ (q 1 - 1 5q x+1 )(l-q- x )R n ^(x- W.n. !. q -,.6 \ q 

= (q~ x - jSq x+1 )(l - qi )R n ( n{x) ; a, p, f,S \q ) (19) 
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(q~ x - qa)(l - j6q x+1 )R n \u(x); qa, g _1 /3, 7, qd \ q ^ 

+ (qa - jSq x+1 )(l -q~ x )R n [n{x- l);ga,g- •'./.-../.> I q 

= (q- x - 7 5q x+1 )(l-<jo)n n [ ,,[.r):n. \ ,,). (211) 

where R n is given by (7). 
For the balanced terminating 403 I a < b > c < d 



q;q) (i.e., one of a, b, cord equals 



q n and a ■ b ■ c ■ d — u ■ v ■ w) we also have the following relation: 



(c-u)(l-vc ){wq -1)4^3 



a,b,q 1 c,d 
u, v, q^ 1 w 



q;q 



(, ,/)(] ,,/ 1 )! irq 1 i) f , i ('°' 6 ' C,? 1 ld 

1 u,v,q L w 



q;q 



(c-d)(wq 1 -b)(l-aqw 1 ) 4 3 [ a ' & ' C ' ^ g; g ] 

\ u,v,w J 



(21) 



In terms of q-Racah polynomials, the last relation can be rewritten as 



(q- x - g 7 )(l - (35q x+1 )(a - l)R n L(x + 1); q~ l a, q(3, 7, q~ l 8 \ g) 

+ (97 " 7<5<f +1 )(1 - q-^T 1 )^ ~ l)Rn (v(x); q^a, q/3, 7, q^S \ q^j 

= (g- - 7 <5<f +1 )(« - a/?g" +1 )(l - c^g"")^ f M*)? A 7, 6 I 9) (22) 



(q~ x - qa)(l - pSq x+1 )(j - l)R n \u(x + 1); a, (3, q^j, 6 \ q 

+ (qa - 7^" +1 )(l ~ q- x Hl- X )(l ~ l)Rn U(x);a, (3, g"^, 5 \ g) 
= (q- x - 7 ^ +1 )(7 ~ a(3 q n+1 ){\ - j-'q-^Rn (W «, P, 7, 8 \ q), (23) 



Proof. To prove the first relation for the basic hypergeometric function we ex- 
pand 403 into series in q and perform straightforward computations in every 
term. 
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To obtain (19) and (20) we simply rewrite the relation (18) in terms of 
g-Racah polynomials using their definition (7). 

Next, we observe that gr-Racah polynomials form an orthogonal basis in the 
corresponding l 2 space. Consequently, we can write dual relations for (19) and 
(20) and these are precisely (22) and (23). It is easily seen that two last relations 
are equivalent to just one relation for basic hypergeometric function (21). □ 

Using the last lemma we obtain the following one: 
Lemma 7.4. 



n>0 



wh 



ere 



1 _ K 2„x+N-t 



I _ K 2g2x-t-S+l ' 



1 _ K 2 a x - T+1 

w 1 (x) = q T + N - 1 - t (l-q*- S - N + 1 )-. 1 " ' 



I _ K 2q2x-t~S+l ' 



while Wt t s(%) stands for the weight function corresponding to the parameters 
t,q,K,N,T,S (see Section 4 and Theorem 4-1 for details). 

Proof. First, we substitute the parameters of g-Racah polynomials given by 
formulas (8)-( 11) into the statement of Lemma 7.3. 

We use (19), (20) in cases (8), (9), and we use (22), (23) in cases (10), (11). 

In all 4 cases we rewrite the corresponding relation in terms of orthogonal 
functions f^{x) and get the following 



cU^iv) - \r¥^Tv^y - ^ - v + \ ^^^y)^y) ( 24 ) 

Multiply the last relation by f„(x) and sum over all meaningful n. 
Since functions /^(y) form an orthonormal basis in the corresponding l 2 
space, 

n 

and the needed relation follows. □ 

Proposition 5.2 implies that the transition probabilities P t +(X,Y) have a 
determinantal form. The last lemma yields that this form is exactly the one 
required for the application of Proposition 7.1. Thus, the theorem is proved. □ 

Applying Proposition 7.1 for the process X(t) we obtain the following state- 
ment. 
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Theorem 7.5. 

Pn {pl ; ^1; • • • ; 7 *^n) — dct [-/^(/c^ , , A^^" , Xj )] ^ n 7 

where 

N-l 

K(k,x;l,y)=J2 -ij;fi(x)f!(y),k>k 

i=0 C i 

K(k,x;l,y) = -J2 c i' l fi(x)fi(y), * < i; 

c*' fc = 1, of = c* • ^ c^" 1 

and functions fi(x) and numbers c\ are given by the formulas (16) and (17). 



7.2 Inverse Kasteleyn matrix 



Let us present another way to view the correlation kernel derived in the previous 
section. 

Recall that we deal with lozenge tilings of a hexagon. Divide every lozenge 
into two unit triangles and color the resulting triangles into black and white 
(west triangle is black) . In this way a tiling turns into a perfect matching of the 
part of the dual hexagonal lattice that fits in our hexagon. Correlation functions 
of the perfect matchings can be computed using Kasteleyn's theorem (see [Ka]). 
Let us describe it. 

Associate to every triangle the midpoint of its vertical side. Note that in 
this way both black and white triangles can be parameterized by the points of 
the two-dimensional lattice. Thus, we can use our usual coordinates (t,x) for 
the triangles. 



OR 




The Kasteleyn matrix Kast(i, x; r, y) is a weighted adjacency matrix. Here 
(t, x) stand for the coordinates of a white triangle and (r, y) stand for the coor- 
dinates of a black triangle. In our case, 



Kast(i, x; r, y) 

Set 
g(t,x) 



nq 

1 
1 





-S/2+x-t/2+l/2 _ 



1 



-S/2+x-t/2+l/2 : 



nq 



(t,x) = (r,y), 

(t,x) = (r-l,y-l), 
(t,x) = (r- l,y), 
otherwise. 



f]\t+x K -tqX(T+N-t-l)+t(S/2-l/2)+t(t+l)/4fj _ K 2q2x-t-S+l-j 



y/w t ,s{x) (q l ;q 1 )s+N-i-x(q;q)T-s+x-t(n 2 q x T+1 ;q) T +N-t 



(25) 



3.3 



Theorem 7.6. Consider n lozenges enumerated by pairs of triangles ((ti,Xi), (r^yi 
The probability that a random tiling contains these lozenges equals 

n 

JjKast(ti, a*; rj,s/i) • det [if ea;t (r l ,y i ;^,a; J )] . j=:L ^, 

i=l 

where 



](t,x) 

the function g is given by (25), and K(r,y;t,x) is given in Theorem 7.5. 
Proof. Kasteleyn's theorem states that the probability to find lozenges 

((ii, xi); (r*i, yi )),... , ((i n , a; 
can be expressed via the inverse of the Kasteleyn matrix: 

n 

JjKast(£i, Xi\ri,yi) -det [Kast _1 (r i; yf, tj, xj)] ij=1 ^ ^ 

4 = 1 

We can compare this statement with Theorem 7.5. Note that Theorem 7.5 
describes the correlation functions of the particles. Consequently the correlation 
kernel 

K(t, x; r, y) = Sfyjfi - K(t, x; r, y) 

(where K (t, x; r, y) is the correlation kernel of Theorem 7.5) is the correlation 
kernel of holes or, equivalently of horizontal lozenges. (See, e.g., [BOO, Ap- 
pendix A. 3] for some details on the particle-hole involution.) 

Since the correlation kernel appears only in a determinant, it is only deter- 
mined up to conjugation: the transformation 

K(t, x; r, y) i-> 9 ^ ,X K(t, x; r, y) 
9{r, y) 

does not change correlation functions. We conclude that 

K(t,x;r, y) 



K q-S/2+x-t/2+l/2 _ ( K q-S/2+x-t/2+l/2-j-l 

should be (perhaps, after some conjugation) the inverse Kasteleyn matrix. 
Let us verify this fact and find the appropriate conjugation factor. 
We have 

\ - 9(t,x) K(t,x;h,z) _ Ar,y) 

gjh, z) Kq -s/2+x-t/2+i/2 _ ( Kq -s/2+ x -t/2+i/2yi ^ as H^ z , r, y) o [t x) . 

(26) 
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First, suppose that t < r — 1. In this case the ((t,x); (r,y)) matrix element 
of the right-hand side is zero while the one of the left-hand side of (26) is 

(tl X ) rtf \ t.l — 1 



K q-S/2+t-x/2+l/2 _ t R q-S/2+t-x/2+l/2\ 



i>N 



( -S/2+2/-r/2+l/2 1 A r-1 fr(\ , r - X ,, , r _i. 



g{r,y) g(r-l,y) g(r-l,y-l) 

Let us find such function g that for every i 

-S/2+y-r/2+l/2 1 \ r-1 fr(\ , r -l,s j-r-Xt i\ 

1 Kg-s/i+v-^+iy ) % JAV) fr fl 1 (y-l) = Q 

9(r-l,y) g(r-l,y-l) 
(27) 



We know that (see Lemma 7.4) 



c*/* +1 ( y ) - \M^- W i(y-i)fKy-i) - J^L^vMy) = o, 
V wt+i.sU/j V w «+i,s(y) 



(28) 



where 



1 _ K 2 x+N-t 



2 n 2x-t-S+l ' 



1 — K 2 q 

i ^.2 z-T+l 

^(x) = ^"^(l - g -(s+"-i-*) ): 1 - '■ 



I _ K 2g2x-t-S+l ' 

while w tt s( x ) stands for the weight function corresponding to the parameters 
t, q, k, N, T, S (see Section 4 and Theorem 4.1 for details). 

For every pair (r, y) we get the following three equations defining g 



9(r,y) oc 



-S/2+y-r/2+l/2 1 

„„-S/2+«-r/2+l/2 



g(r - 1,3/ — 1) oc - 
g(r- l,y) oc 



V w r,s(y) 

1 



\/wr-i,s(y - 1 (y-l)' 

1 



where the proportionality coefficient is the same for all three equations (but it 
may depend on the pair (r,y)). One checks that g given by (25) satisfies these 
relations and after the conjugation with g the ((t, x); (r, y)) matrix element of 
the left-hand side of (26) is zero. 

Next, suppose that t > r. In this case the ((t,x); (r,y)) matrix element of 
the left-hand side of (2 ) is zero by the similar reasoning. 

If cither t — r or t = r — 1 the argument becomes a little more involved, but 
the computation requires no new ideas. 

□ 
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8 Bulk limits. Limit shapes. 



8.1 Bulk limit theorem 

In this section we compute so-called "bulk limits" of the correlation functions 
introduced in the previous section. 

We are interested in the following limit regime. Fix positive numbers S, T, 
N, t, x, q. Introduce a small parameter e«l, and set 

S' = Se- 1 + o(e- 1 ), T = Te- 1 + o(e~ 1 ), N = Ne" 1 + o^ 1 ), q = q e+ °^ . 

Consider also integer valued functions ti = U(s) and Xi = Xi(e), i = 1, . . . , n, 
such that 

lim eti(e) = t, lim exi(e) = x, i = 1, . . . ,n, 

e^O e— »0 

and pairwise differences U — tj, and Xi — Xj do not depend on s. 

Then the correlation functions p n computed in Theorem 7.5 tend to a limit 
p n which depends on the parameters of the limit regime q, S, T, N, t, x and the 
original parameter k. 

We consider the region where the limit correlation functions are nontrivial. 
This region is commonly referred to as the "bulk" , sometimes also called the 
"liquid region" . This is a simply connected domain inside the hexagon. 

The main result of this section is Theorem 8.1. 

Note that the first limit correlation function allows us to predict the limit 
shape which appears in our model. 

Theorem 8.1. We have 



lim p n (h, xr, ■ ■ ■ ;t n ,x n ) = det K(ti,Xi\tj,x 3 ) 



i,j=l,...,n 



wht 



K(x, s;y,t) = — <j) (1 + cw) t ~ s w^^dw. 



Here the integration is to be performed over the right side of the unit circle when 
s > t and over the left side otherwise, 



q T - 2t (l - q-( 5+N - x ))(l - q x ) (1 - At 2 q x+N - S )(l - K 2 q x - T ) \ 5 
(1 - q x + T - t - s )(l - q-t- N + x ) (1 - K 2 q x + N - t )(l - ^V=i=S) J 

and (j) is given by the formula: 

q- N (l - q N )(l - q- J - N )(l - K 2 q- t - s + 2 «) 2 + A + B 

a> = arccos ; , 

2VAB 

where 

A = (1 - q- S - N+x )(l - ^ 2 q- T+x )(l - q^^l - ^q^ 5 ^), 
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B = q- 2N - T (l - q x )(l - K V t+N + x )(l - q- t - S+T+x )(l - K 2 q- S+N+X ). 

// the expression under arccos is greater than 1, then we set <p = 0. If the 
expression is less than —I, then <f> — tt. 

Setting v = —cw (and omitting some "conjugation factors" again) we get 
the incomplete beta-kernel form of the integral, cf. [OR] , 

1 

K(x,s;y,t)= — <b (1 - u)*" s v^^dv. 



c-e~ 



Here the contour of integration intersects (— oo, 0), if s > t, and intersects (0, 1) 
otherwise. For an explanation of the relation of the incomplete beta-kernel and 
Gibbs measures see [KOS], [BS]. 

It is not hard to compute that z = —ce %< ^ has the form 



1 q T + N - 1 

z 



2(1- q x + T - t - s )(l - q- t - N + x )(l - K 2 q- t + N +x)(l - «;2q*-t-S) 
q- N (l - q N )(l - q- T - N )(l - JcT*-^ 3 *) 2 + A + B 



i\/4AB - (q- N (l - q N )(l - q- J - N )(l - K 2 q -t-S+2x)2 + A + B) 2 



with 

A = (1 - q- S - N+x )(l - K V T+X )(1 - q- t - N+x )(l - ^q-*" 5 ^) 

and 

B = q- 2N - T (l - q x )(l - K 2 q- t+N+x )(l - q- t - S+T+x )(l - K 2 q- S+N+X ). 

Proposition 8.2. The parameter z defined above coincides with the one defined 
in Theorem 2.1. 

Proof. The quadratic equation satisfied by z is 

Z 2 - (z + z)Z + zz = 0. 

Substituting the expression for z given above, one obtains a relation equivalent 
to the one in Theorem 2.1. □ 

8.2 Proof of the bulk limit theorem 

In this section we prove Theorem 8.1. 
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Recall that the correlation kernel before the limit is given by 

N-l 

K(x, k;y,l) = -J2 Ci' l fi(x)fi(y), * < h 

i>N 

^ = 1,^ = ^-^ c\~\ 

The functions f^(x) and the coefficients were defined in Section 7. 

First, let us consider the case k = I. We want to find a limit of the projection 
kernel 

JV-l 

V t (x,y)=J2fn(x)fn(y)- 

71=0 

In order to find a limit of Vt{x,y) we use the spectral projection method 
proposed by G. Olshanski, see [B02] and [O]. 

We want to consider V^x^y) as a matrix element of the operator Vt- It 
turns out that finding the limit of the operator is easier than computing the 
limit of the matrix elements. Note that functions f^ix) are eigenvectors of 
some difference operator (it will be explicitly given below). The projection 
operator can be regarded as the spectral projection on the segment containing 
the first N eigenvalues of this difference operator. Now, to find the limit of the 
spectral projection operators we will take the limit of the difference operators. 
Note that both the difference operator and the spectral segment are varying 
simultaneously. 

To justify the limit transition we use some facts from functional analysis. 

Consider the set Z°(Z) of the finite vectors from ^(Z) (i.e., the algebraic span 
of the basis elements S x ) as a common essential domain of all considered differ- 
ence operators. It will be clear from the following that the difference operators 
strongly converge on this domain. It follows that the operators converge in the 
strong resolvent sense (see [RS], Theorem VIII. 25). The last fact, continuity 
of the spectrum of the limit operator, and Theorem VIII. 24 from [RS] imply 
that the spectral projections associated with the difference operators strongly 
converge on the set of finite vectors to the limit spectral projection associated 
with the limit difference operator. 

Now we present some details and computations. 

Note that since q-Racah polynomials are eigenfunctions of a certain differ- 
ence operator (see [KS]), the same is true for the functions fn{x). The difference 
operator is 
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<T"(1 - - a/3<f +1 )/^) = B(x)f n (x + 1) 1 



' w tiS (x) 
wt,s(x + 1) 



- [B(x) + D(x)]f n (x) + D{x)f n {x - 1)J m f {x) (29) 

where Wt t s(x) is the weight function corresponding to the parameters t, q, k, N, T, S 
(see Theorem 4.1), and 

(1 - a q x+1 ){l - (36g x+1 )(l - lq x+1 ){l ~ l5q x+1 ) 
B(x) = 



D{x) 



(1 -7^ 2a; + 1 )(l -j6q 2x + 2 ) 

qQ. - g x ){a - <y6<f)(P - jf)Q. - Sg x ) 
(1 --fSq 2x )(l - 1 5q 2x + 1 ) 



(Here a, (3, 7, 8 are the corresponding parameters of g-Racah polynomials.) 
We find 

q 2N+T-U 1 _ K 2 q 2 x -t-S+Z ) 

wt, s (x + l)M, S (or) - - _ ^ 2g2rc _ t _ s+1 

(1 - q x - t - N+1 ){l - q x - S - N+1 ){\ - K 2 q x - T+1 ){l - it^'- 5 + 1 ) 
X (1 - q x+1 ){l - g T-S-t+x+l)( 1 _ K 2 ga; +JV- t +l)( 1 _ ,,2^+JV-S+l) ' 



w tj s(x - l)/w t ,s(x) 



1 _ K 2^2a;-t-S+l 

(1 - g a: )(i - g T - s - t+ ^)(i - kV+^-'xi - ^V +jV - s ) 

X (1 - - ^- S - Ar )(l - K 2 q*- T )(l - ^qx-tS) ' 

Substituting the g-Racah parameters of our model (see Section 4) into B(x) 
and D{x) one computes 



B(x) 



D(x) 



(1 - qS-N+x+l^ _ K 2 g -T +3: +l )(1 _ g-t-N+x+1^ _ ^-t-S+x+lj 
(1 _ K 2g-t-S+2x+l^^ _ K 2g-t-S+2a;+2^ ' 

_ q(l - g^(g-^ - ^V^+'Xg 5 - 7 ^ - g^^!Kj - K 2 q- S+N + X ) 
~ (1 - K 2 q -t-S+2x)(i _ ^-t-S+22+1) 

.1-2JV-T (1 - «*)(! - ^ 2 ^ t+A,+X )(l - g-*- S + T +^)(l - K 2 q- S + N+X ) 
(I- n 2 q- t - s + 2x )(\- K 2 q- t - s + 2x + 1 ) ' 
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X 



It follows that 
B(x) 

^w t ,s{x + l)/w t ,s{x) 

_ ^ q l-2N-T {1 _ qX +l) {1 _ q TS-t+ x+ l ){1 _ q S-N +x +l ){1 _ q -t-N +x +l ) 

(1 _ _ _ K 2^+JV-t+l )(1 _ K 2 q x+N-S+1} 

(1 — K 2 q~ t ~ S+2x+1 )(l — K 2 q 2x ~ t ~ S+3 )(l — K 2g-i-S+2x+2^2 

D(x) 



\fw t .s{x - l)/w t ,s(x) 

= {q l - 2N - T {\ - <f )(1 - q x - t - N )(l - q*- S - N )(l - ^- t ~ s + T )) i 
'(1 - K 2 q- t+N + x )(l - K 2 q- S + N+X )(l - K 2 q x - T )(l - n 2 q x -^ S ) 



(1 - n 2 q- t - s + 2x + l ){\ - kV* - * -5-1 )^ - n 2 q- t - s + 2x ) 2 
The eigenvalues in the left-hand side of (29) become 

q- n (l ~ q n ){l - q -T-2N+n+i^ n = , . . . , JV - 1 

Note that in the "bulk limit" regime, q — > 1, N, T, S — > oo as e — > in such 
a way that q N , q T , q s have finite limits. 

In the limit , B ^ tends to some constant and , SMI tends to 

\/ w(x+l)/w(x) \l w(x— l)/w(x) 

the very same constant. After dividing by twice this constant the limit operator 
becomes 

/M „ /(, + l) + /(,-l) _ /(i) M + B 



2VAB 
where 

A = (1 - q- s - N+x )(l - K 2 q- T+X )(l - q- t " N + x )(l - K 2 q- t_S+X ) 

and 

B = q- 2N - T (l - q x )(l - K 2 q- t+N+x )(l - q- t - S+T+x )(l - K 2 q- S+N+X ) 
while the spectral interval becomes 

V N (1 - q N )(l - q- T - N )(l - «V t-S+2 *) 2 



2VAB 



,0 



Since the spectrum of the operator / — > £ig±ll±Z(£ 1 ) j s while 
2v^fg ^ ^' one can ec l u i va l en tly write the operator in the form 

> /(s + l) + /(a:-l) 
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with spectral interval 

NV1 _ n-T-N^-l _ *2 n -t-S+2x^2 _|_ 4 _|_ R 1 

(30) 



q- N (l - q N )(l - q- J - N )(l - K 2 q ~ t-S+2x ) 2 + A + B { 



2VAB 

If we perform the Fourier transform ^(Z) — * L2{S 1 ), where S 1 is the unit 
circle in C, we get an operator of the projection on the part of the unit circle 
with .T-coordinate varying over precisely the spectral interval (30). 

If we do the inverse Fourier transform we will get the discrete sine kernel, 
cf. the end of Section 3.2 in [Gor]. 

Next, let us consider the case k < I. The prelimit correlation kernel is given 

by 

K(x, k;y,l) = -J2 'fi(x)fi(y), * < h 

i>N 

c*' k = l, a? = c\~\ 

Let us decompose the correlation kernel (i.e., the operator given by it) into 
the product of the static projection kernel: 

N-l 



i>N 

and a collection of transition operators (or their inverses) Uh(x,y) with 

f E 4ft(x)^ +1 (y), xeX hl ye X h+1 , 
U h (x,y) = <^ i>o 

I for other x, y. 

For the operators P't(x, y) we can use the same methods as for Vt(x, y). We 
get minus the operator of projection on the part of the unit circle complementary 
to the spectral interval (30). 

Let us turn to the transition operators Ut- 

By virtue of already proved facts (see Section 7) , we obtain 

U t (x,y) = const- / w t,s(x) \ Wl ^gv +w (x 
y Wt+i,s\V) 

where 



x(2N+T-l)M _ K 2 2x-t-S+l\ 



(q;q)x(q;q)T~s~t+x(q 1 ;q x ) t +N- x -i{q u ,q 1 )s+n-x-i 

i 



(kV- t+1 ; «)r+Ar-t(«V-*- s+1 ; q)N +t ' 
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1 _ K 2 x+N-t 
W (x) = -(1 - 1 — q 



I _ K 2g2x-t-S+l ' 



and 

Thus, 
where 

w (x) = 
and 

wi(x) = 



2„x-T+l 



U t (x,y) = const - [wi(x)S% +1 + w (x)5%] , 



(1 - g 3: + T - t - s )(l - q- t - N+x )(l - K 2 g "- t - s )(l - K 2 q x+N - t ) 

(1 - K 2 g 2 a; -t-S+l)( 1 _ K 2 q 2x-t-S) 



(1 - ^^-t-ST+l)^ _ K 2 9 2x-t-S+2) 

Passing to the limit we get the operator 

U(ar,») = consi- [U x ^ +1 + U S«] , 

where 

(1 - q x+T - t - s )(l - q- t - N+x )(l - rcV~ t_S )(l - K 2 q x+N ^) 



U = 
and 

Ui = 



(1 - /c2q»c-t-S) (1 _ K 2 q 2x-t-S) 
q T - 2 *(l - q x - s - N )(l - K 2 q x - T )(l - q x )(l - K 2 q x+N - S ) 



(1 - /t 2 q 2x - t - s )(l - « 2 qa»-t-S) 

Equivalently, 

U(af, y) = consi • [u<5* +1 + 8»] , 

where 

_ Ui _ / q T - 2t (l-q x - s - N )(l-q x ) (1 - K 2 q x+N - S )(l - K 2 q x - T ) 
U ~ ~ ^(1 - q x + T - t - s )(l - q-t- N + x ) (1 - K 2 q x + N - t )(l - «;2qx-t-S) 

Fourier transform gives us the operator of multiplication by const(l + u/w), 
where w is the coordinate on the circle \w\ = 1. 

If we now multiply all necessary operators, perform inverse Fourier transform 
and substitute w — > 1/w in the resulting integral, we get the desired limit kernel. 
Note that the constant prefactor in U can be omitted since it corresponds to 
the conjugation of the kernel that does not affect the correlation functions. 

The final case k > I is similar. The interested reader can find some details in 
jGor] where similar computations (with Hahn polynomials instead of g-Racah 
polynomials) were performed. 
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9 Computer simulations. Different limit regimes. 



Using the perfect sampling algorithm described in Section 6, we performed some 
computer simulations; the program that we used can be found at http://www. 
math.caltech.edu/papers/Borodin-Gorin-Rains.exe. We are mostly inter- 
ested in the case when the hexagon is large, since in this case we can see some 
limit shapes appearing. In all pictures we color the three types of lozenges in 
three different colors (as in Figure 1). When we draw big pictures, we erase the 
borders between lozenges and get some coloring of a hexagon in three colors 
which can also be viewed as a stepped surface in K 3 . 

Although we mostly show not very big pictures, our algorithm can generate 
random tilings of a 1000 x 1000 x 1000 hexagon in a reasonable amount of time. 

The following picture shows a plane partition in a 70 x 90 x 70 box sampled 
from the distributions with parameters q = 0.97, k — 1. The formation of a 
limit shape with frozen regions is clearly visible on the picture, and the next 
picture shows the border of the frozen region as predicted by Theorem 2.1. 




By the appropriate limit transition, we get plane partitions distributed as 
^volume _ f rj owm g picture shows a random plane partitions in a 70 x 90 x 70 
box and the corresponding theoretical frozen boundary with q = 1.04. 
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As was explained in Section 4, if we send q — > 1 in the original model, 
then the weight of a horizontal lozenge becomes a linear function in the vertical 
coordinate. If we tune the parameters in such a way that our linear function has 
a zero at the bottommost point of a hexagon, then we get the following random 
plane partition in a 70 x 90 x 70 box, cf. the corresponding theoretical frozen 
boundary. 




We see that the border of frozen region has a node near the bottommost 
point of a hexagon. 

Another interesting case to consider is when q does not tend to 1 as the 
size of the hexagon tends to infinity. Look at the following picture where the 
random plane partition with parameters q = 0.9, k = 1 in a 70 x 90 x 70 box is 
shown. We see that the surface becomes different from the ones shown above. 
We call this new class of surfaces waterfalls. 



44 



The next picture shows an even more degenerate case q — 0.7. 



We hope to study the asymptotic behavior of these cases in a later publication. 

Finally, let us turn to the trigonometric g-Racah case. In this case the weight 
of a horizontal lozenge is sin(a(j — (5+l)/2) + /3). One can tune the parameters 
a and j3 in such a way that this weight becomes zero at both the topmost and 
the bottommost points of the hexagon. The boundary of the frozen region has 
two nodes in this case. 
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10 Appendix. Lozenge tilings and biorthogonal 
functions 



We were led to consider the g-Racah and Racah cases via the realization that, 
in much the same way that uniform lozenge tilings are related to Hahn polyno- 
mials, so were g-weighted lozenge tilings related to g-Hahn polynomials. This 
naturally led to the question of whether more general (discrete) hypergeometric 
orthogonal polynomials could arise in this way. In fact, as we will consider in 
this section, one can generalize even further, to the elliptic analogue, certain 
biorthogonal elliptic functions due to Spridonov and Zhedanov [SZ]. 

One way to derive the required lozenge weights (essentially equivalent weights 
were considered by Schlosser in [Sch], although the connection to plane parti- 
tions or lozenge tilings was not made explicit there) is via the following desider- 
ata: 

First, we want the total weight of all tilings of a hexagon to be "nice", in 
the sense that it should be expressible as a product of simple theta functions 

e p (x) :=n(i-A)(i-p**7*); 

the same should apply to the individual weights as well. Note that, as is tradi- 
tional in much recent work on elliptic special functions, we use the multiplica- 
tive form for our elliptic curves and theta functions. This can be translated to 
the usual doubly-periodic form by composing with the singly-periodic function 
x i— > exp(27nx), but the multiplicative form makes certain degenerations simpler 
to obtain. 

Next, the sums that arise should be hypergeometric, in the sense that any 
parameters should vary along geometric progressions as one moves around in 
the hexagon. 

Finally, the form of the weight of a given cube (i.e., the ratio of the weights 
of the two ways to tile any given unit hexagon) should be invariant under the 
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symmetry group of the triangular lattice, except that certain reflections should 
invert the weight. Note that the weights of cubes are gauge-invariant, and any 
set of choices of cube weights corresponds to a unique gauge-equivalence class. 

We may rewrite these criteria in terms of plane partitions, noting that the 
requirement that the cube weights come from lozenge weights places the addi- 
tional restriction that if w{x, y, z) is the weight of a cube with centroid (x, y, z), 
then 

w(x, y, z) = w(x + l,y+l,z+l). 
If we consider plane partitions inside alxlxn box, we thus require that 



0<i<n 0<m<Z 



II «>(l/2,l/2,m + l/2) 



should be nice, and should indeed correspond to a hypergeometric sum. The 
isotropy condition then implies more generally that 

IT w{x,y,z + m) 

0<l<nO<m<l 

is a hypergeometric sum for any half-integer vector (x,y,z). In particular, this 
sum is doubly-telescoping, in the sense that the sum over any subinterval is 
nice. In particular, if we assume that the sum should be a special case of the 
Frenkel-Turaev sum [FT], we find that it should have the form (up to some 
convenient changes of parameters) 

q'epjabq 21 - 1 ) O p (a/q,a,b/q,b) = fl p ( g "- 1 a6,g"+ 1 ,a,6) 

where a and b depend on (x, y, z). In particular, we find 

qe p (q m - 1 a,q m - 1 b,abq 2m+1 ) _ q%{q m - 1 a,q m - 1 b,q- 2m - 1 /ab) 



w(x, y, z+m) = 



9 p (q m+1 a,q m+1 b,abq 2m ~ 1 ) 6 p (q m+1 a,q m + 1 b,q 1 - 2m /ab) ' 



where a, b depend on (x,y,z); consistency then implies that q~ z a and q~ z b are 
independent of z. 

Rotating the picture by 120 degrees gives a similar expression for w(i,j, k) 
with the dependence on i or j factored out; comparing the results leads us to 
an expression of the form 

_ q z p {qy +z - 2x - 1 u 1 ,q x+z - 2 y- 1 U2 1 q x+y - 2z - 1 U3) 

W \ X >y-> Z ) - p ( q V+ Z -2x+l Ulj q x+z-2y+l U2j q x+y~2z+l U3 ) > 

where u% = a, 1*2 = b, U3 = 1/ab; i.e., u\, U2, U3 are generic such that U1M2M3 = 
1. We will see that this indeed gives rise to a factored sum over plane partitions 
in a cube. Note that if we rewrite this expression in terms of new variables 
til = q v+z ~ 2x Ui, u 2 = q x+z ~ 2y u 2 , u 3 = q x+v ~ 2z u 3 , then 

q 3 9 p (u 1 /q,u 2 /q,U3/q) 

w(x,y,z) = — f— — — , 

p {qu 1 ,qu 2 ,qu 3 ) 
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which can be described in a canonical way: w is the value at q of the unique 
elliptic function with simple zeros at Hi , simple poles at u~ 1 , and the value 1 at 
1. It follows, in particular, that w is invariant under shifting any parameter by 
p, as well as under all modular transformations. 

To write this in terms of lozenge weights, there are, of course, a number 
of choices one could make. One convenient choice is to allow only horizontal 
lozenges to have nonzero weights. We must thus have 

w(i,j + l) _ qe p {q^- M / 2 - 1 u l ,q3+ :il ' 2 - 1 U2 1 q 2j + 1 u l u 2 ) 
~~ 9 p (qj- 3i / 2 + 1 u 1 ,qj+ 3i / 2+1 u 2 ,q 2 i- 1 u 1 u 2 ) 

where w(i,j) denotes the weight of a horizontal lozenge with upper corner at 
(recall that in terms of 3-D coordinates, i = x — y, j = z — (x + y)/2). 
This recurrence is straightforward to solve, and one obtains 

«,(U)= C (i) - q j ~ 1/2 (uiu 2 ) 1/2 e p (q^- 1 u 1 u 2 ) 



e p (qi-^/ 2 - 1 u 1 ,q3-^/ 2 u 1 ,q3+^/ 2 - 1 u 2 ,qi+^/ 2 u 2 ) ' 

where C(i) is an arbitrary (non- vanishing) function of i. 

If we view (with a mind to applying Kastclcyn's theorem) the lozenge weight 
as a matrix indexed by a right-pointing and a left-pointing triangle, we find that, 
coordinatizing triangles by their upper corners, 

»((«), («)) = CW (»«>"V-^<^„ 2 ) 



Jp(qi- 3t / 2 - 1 u l7 q3- 3i / 2 u 1 , q i+ :ii / 2 - 1 u 2 , q i+ 3i / 2 u 2 ) 
w((i,j), (i + 1, j + 1/2)) = 1, w((i,j), (i + l,j - 1/2)) = 1, 

and all other values are 0. 

Let H%°'^\ represent the parallelogram 

x <i< a;i,2/o < j + V 2 < 2/ii 

and observe that the restriction of w to triangles in is a square matrix, 

and we can thus attempt to invert w in such a region. In fact, not only can we 
explicitly invert w in such parallelograms, but the result is independent of the 
choice of parallelogram. 

Theorem 10.1. The inverse transpose of w in has the form 

W((i Q ,jo),(ii,ji)) = S io<il [] C{k)-{ Ul u 2 )^-^l 2 

i <fc<ii 

v A . . , n ^. . ,„f_lVo-io/2-ji+ii/2-l„(ii-io-l)(ii-io+4ii-2)/4 

x o n +i 1 / 2 < 3a +i„/ 2 [ ij q 

fl p (gJo+W2-Ji -«i/2+l, q 3Q+*o/2+n^i/2 UlU2; g) jl _ jg _i 

9 p (q, q^- 10 / 2 -^^, (7-J'i- 34 i/ 2+1 Mi, gJ'i+^+'i/ 2 ^, ^'o+3»o/2+l U2; q) ll _ 4o _ 1 ' 

w/iere 

9 p (x;q) k ■= }[ d P (q l x). 

0<i<k 
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Proof. Note that by Kastelcyn's theorem, W is the total weight (up to sign) 
of all lozenge tilings of the parallelogram that omit the two given triangles. In 
particular, W((io, jo), (*i, ji)) must vanish if i > i\ or j + io/2 < ji + ii/2 
simply because then no such tiling exists. 

Now, the claim that W is the inverse transpose of w reduces to the statement 

W((to, jo), (ii, ji)M*i, ji) + W((io,jo), (*i + 1, ji - 1/2)) 

+ W{{i ,j ), (*i + 1, ji + 1/2)) = fyojoMiiJi)- 

This holds trivially for «o > «i, since all three terms vanish, and similarly for 
*o = h, ji > jo- When i = ii, j = j 1; the claim reduces to 

W((io,io),(io + l,io-l/2)) = l, 
while for «o — Jo > ji, w ^ have 

W((to, j ), (io + 1, ji - 1/2)) + W((i ,io), (io + 1, ji + 1/2)) 

= (-l)J'o-Ji-l + (_!)Jo-ji =Qj 

as required. 

It remains to consider the case zo < ii- If ji+h/2 > jo + io/2, then all three 
terms again vanish, while if ji + i 1 /2 — j + io/2, the third term vanishes, and 

W{{i ,j ), (*!, h))w{h, ji) + W((i , jo), (*i + 1, ji - 1/2)) = 

as required. Finally, when ji + ii/2 < jo + *o/2, so that all three terms survive, 
if we divide by 

W((i , j ), (ii, ji))w(ii, ji) 

p (q i i- i o,qio- l o/ 2 - i i- 1 u 1 ,q^>+' l o/ 2 +^U2,q 2 ^- 1 u 1 u 2 ) ' 

we simply obtain a special case of the addition law for 9 p , in the form 

9 p {a z 1 a\z, a 2 z, a a 1 a 2 /z) — 9 p {a a a ll a a 2 , aia 2 , z 2 ) 

+ 6 p (z/ao, ai/z, a 2 /z, a a 1 a 2 z)za n = 

with 

ao = g -J'o/2-3Jo/4+ji/2+3ii/4^ ^ _ ? jo/2-io/4+ji/2-5ii/4-l Wi 
02 = q j ^ 2+3i "/ 4+: > 1 / 2+3il ^u 2 , Z = flJo/2-io/4-ji/2+i 1 /4_ 

□ 

Remark. One major source of guidance regarding the form of W is that it 
corresponds to an enumeration of plane partitions in a rectangular parallelepiped 
with dimensions m x n x 1, say; i.e., a sum over ordinary partitions. If we first 
sum over the first part of the partition, we find that W should look like the term 
of a (singly) telescoping hypergeometric sum. One can also, of course, compute 
"small" values of W in the case p = 0, and look for patterns in the resulting 
factorizations. 
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Lemma 10.2. Consider the domain 




i ( i 0+ l,j-I/2-x k )f ki 



i +H-I/2-b-c) 



The total weight of lozenge tilings of this domain is equal to a constant 
independent of {xk} times 

yr (-l) x %(q I+1 , g J -^+ 3 '°/ 2 M, g- J + 2 -JQ-3Wy M2 , q'+^/u^q)^ 

i<i< c e p( q ' q 2I ~ jo+Mo/2 /ui,q 2 - jo - 3i ^ 2 /u2,q 1 - 2 ^/u 1 u 2 ;q) xl 



times 



Ul<k<Kc q- Xk P (q Xk - Xl , q Xk+Xl+I - 2la + 1 /u lU2 ) 



rii<;<c p (q 1 - I +3«- 3l «/ 2 -^u 1 ,q x '-^- :ii o/^/u2:q) c -i ' 
Proof. The problem is equivalent to computing the weight of tilings of the do- 



main 




(i +U -I/2) 



(i +l,,jg-I/2-b-c) 
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According to Kasteleyn's theorem, cf. [Ka], we obtain 

det [W((i - k,j - k/2 + 1), (to + I, Jo - 1/2 - ari))]*,j=i • 

Now, we can write 

W((i - kjp - k/2 + 1), fa + 7, jo - 1/2 - Xl )) 
W((i - 1, jo + 1/2), fa + /, jo - 1/2)) 

= W{{ip - 1, jo + 1/2), (tg + 7, jo - 1/2 - a;,)) 
W((io - 1, jo + 1/2), fa + I, Jo - 1/2)) 

W((i - k,j - k/2 + 1), fa + I, jo - 7/2 - sQ) 
X W((to-l,j + l/2),(t +/,io--f/2-a;j)) ' 

noting that (for generic values of the parameters) W((io — 1, jo + 1/2), (io 
I, j -1/2)) ^0. Since 

W((ip - 1, j + 1/2), fa + I, jo -I/2- Xl -l)) 
W((i - 1, j + 1/2), fa + 7, j - 1/2 - xi)) 
e p (q x 'q I+1 ,q Xl q I -^+ 3i ^ 2 /ui,q x 'q- I+2 -^- 3ln / 2 /u2,q Xl q I+1 - 2]0 /u 1 u 2 ) 
0p(q x, q,q Xl q 2I ~ JO+3io/2 /u 1 ,q x 'q 2 -i°- 3i o/ 2 /u2,q Xl q 1 - 2jo /u 1 u 2 ) 

we find 

W((i - 1, j + 1/2), fa + I, jo - 1/2 - xi)) 



W((io-l,j + l/2),(i +IJo-I/2)) 
= (-1) 



Xl p (q I+1 , g / - J °+ 3l °/ 2 /ui, q- I+2 -^- 3l «/ 2 /u 2l q'+^o /u x U 2 ; q) xi 



9p(q, q 2I ~~ jo+3i °/ 2 /u 1 , qr2-j n -3io/2 gl-2jo /u\U 2 \q) 3 

For the other factor, we similarly have 

WXfa - k,j - k/2 + 1), fa + 7, jp - 7/2 - Xl )) 
W((i - 1, jo + 1/2), fa + I,jo - 1/2 - xi)) 

6 p {q- x \q l+Xl - 2 ^ +I '/«i« 2 ;g)*-i 



oc 



6l p(9 1_/+io ~ 3io/2 ~ ;r " u i J , ? ;i:! ~ : ' ~ 3lo/2+2 / u 2; £?)fe- 



where we have removed factors independent of , and observed that the right- 
hand side vanishes when < Xi < k as required. As a function of 

= n x '+'/ 2 -^ + l/2 {uiU2) -l/2 



zi = q 



this is invariant under zi i— > and we may thus apply Corollary 5.4 of [Wa] 
to conclude that 



det 



W((»o - k,j - k/2 + 1), (i + I, j - 1/2 - xi)) 



W((i - 1, jo + 1/2), fa + I, jo ~ 1/2 - xi)) 



fej=i 



IIi<fc<i<c q- Xk P (q Xk - Xl , g 35 *+ 3! '+ J - 2 J° +1 /Mm2) 
ril<i<c ^(g 1 - r+io-8io/a -' E, t»l, g»«-J'o-3io/2+2/ U2 . ? ) c _ 1 • 

□ 
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Lemma 10.3. Similarly, the total weight of lozenge tilings of the domain 
(i +I,j -V2).\ 



(i +I,i-I/2-b-c 




(i +a,jo-a/2) 



(i +a,j -a/2-b) 



(i +a-c,j o -a/2-b-c/2) 



is equal to a constant independent of {xk} times 

(-l) x "0 p (q 1 - b - c , g«-J'o+2+3io/2/ Ul) q -a+c-) -Zi /2 J u ^ ^o+a+b+l/^^. ^ 



n 



times 



O p (qI-a-b+l^ q I+2-j +3io/2+a-c/ Ul ^ q -I-j -3io/2 j u ^ ) q I-2j +b+c+l / UlU ^ g) a 



rii<fe<;<c q- Xk 6 p {q x "~ xl . q Xk+xl+I ~ 2jo+1 /u lU2 ) 



]\i <k<c p {q Xk+I - ja+ * la / 2+2+a - c /u l ,q^+ zi «/ 2 + a + 1 - c - x *U2-,q)c~l 
Proof. The problem is equivalent to computing the weight of tiling of the domain 
(i 0+ IJ -I/2)' 

> 

<i +a,j -a/2) 



(i +I,i-I/2-b 




(i +a,j -a/2-b) 



(i g +a-c,j -a/2-b-c/2) 



Again, by the Kasteleyn theorem, we need to compute 

det [W((io + I, jo - 1/2 - x k ), (i + a-c+ l,j - a/2 - b - c/2 + J/2))] . 
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Again, 

W((i + I, j - 1/2 - x k ), (i + a-c+l,j -a/2-b- c/2 + 1/2)) ? 0, 
allowing us to factor the weights accordingly: 

VK(fa + I, j - 1/2 - x k ), {io + a-c+l,j -a/2-b- c/2 + 1/2)) 



(-1) 



W((i + I, jo - 1/2), (i + a-c+l,j -a/2-b-c/2+ 1/2)) 
!h epjq 1 -^, g"-J°+ 2 + 3 '°/ 2 / Ul , g - n + e -J°- 3 '°/ 2 /u 2 , q- 2 3«+ a + b + 1 /u 1 u 2 ; q) Xk 
8 p (q I - a - b +\q I + 2 -3°+ 3l °/ 2 + a - c /u 1 ,q- I -^- 3 *°/ 2 /u 2 ,q I - 2 J°+ b + c + 1 ^ 



and 



W((i + I, j Q - 1/2 - Xkh fa + a - c + I Jq -a/2-b- c/2 + 1/2)) 
W{{io + I, jo - 1/2 - x k ), fa + a-c+l,j -a/2-b- c/2 + 1/2)) 
6 p (q Xk - b - c+ \q 2j °- I - h - c - Xk u 1 u 2 ;q) l - 1 
^ 0p(9 /_JO+3io / 2+2+a_c /ui,(7 : ' o+3io / 2+a+1_c_;l;fc U2;9);-i ' 



Lemma 10.4. The weight of all tilings of the hexagon 
(i„io> 



\ (i n +IJ -I/2) 



□ 




(i +aJ -a/2) 



(i +I,jb-I/2-b-c 



(i +a,j -a/2-b) 
(i +a-c,j g -a/2-b-c/2) 



such that the horizonal tiles along the line i = io+I are as specified, is a constant 
independent of {xk} times 



n 

Kk<c 



e ^ q 2 Xk +I+l-2j j UiU2) q ^e p {q I+1 ,q 1 - b - c ) 



p (q I+1 ^/u lU2 ) 



O^q 1 —^) 



6p{q I - io+ii °l 2 /ui, q- a+c -^- 3i °/ 2 /u 2 , q I+1 ~ 2 ^ /u x u 2 , q- 2 ^ +a+b+1 / Ul u 2 ; q) Xk 
9 P (q I+2 - : > 0+3 ' lo/2+a - c /u 1 ,q 2 -3»- 3i o/ 2 /u 2 , q 1 - 2 ^ /«i« 2 , qt-Vo+b+c+i / Ul u 2 ; q) Xk 
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times 



\<k<l<c 

1 



n 

Kk<c 



9 p (q^k+I-J +3 l0 /2+2+a-c / u ^ q I-j„+3t„/2+2+a-c / Ul ;q) c -1 



n 



9 p (qjo+3i /2+a+l-c-x kU2 ^ qja+3i /2+a+l-c-x k u ^ . g) c _ 1 



Remark. The first "univariate" factor is just the product over the variables of 
a special case of the weight function considered by Spiridonov and Zhedanov 
in [SZ], and the two sets of elliptic hypergeometric biorthogonal functions con- 
structed there are triangular with respect to the functions appearing in the two 
determinants. We could thus replace the second "cross-term" factor by a prod- 
uct of Vandermonde-style determinants in which the fc-th row consists of the 
(fc — l)st biorthogonal function evaluated at x\ through x c , analogously to our 
calculations for the q-Racah case above. 

Using either lemma, we can in particular obtain a formula for the total weight 
of all tilings of a hexagon. Written in terms of plane partitions, we obtain the 
following elliptic analogue of MacMahon's identity. 

Theorem 10.5. Let p, q, u\, u 2 , 113 be generic parameters such that \p\ < 1, 
U1M2U3 = 1. Then 

^ 6Jqi+ k - 2i + 1 Ui,q i + k - 2 i+ 1 u 2 ,q i +i- 2k + 1 u 3 ) 
ncaxbxc(i,j,k)en pvy ' y ' y J/ 



q abc Yl » P W +]+k -\q 3+k - l - 1 u l ^ +k -l- 1 u 2l q^- k - 1 U3 } 



9 p {q*+3+ k - 2 ,q3+ k ->u liq *+ k -lu 2l q*+l- k U3) 

Proof. The term corresponding to a specific plane partition II is easily seen by 
induction to be the ratio of the weight of the corresponding tiling to the weight 
of the tiling associated to the empty plane partition. The claim follows by 
simplifying the corresponding determinant of W using Corollary 5.4 of [Wa]. □ 



Although it is possible to arrange for the elliptic weights to be positive, there 
are difficulties in the analysis. For one thing, the algebra required to replace 
orthogonal polynomials by biorthogonal functions in constructing the kernel 
has not been fully developed. A further complication in computing the limit 
kernel is that, although the biorthogonal functions do satisfy reasonably simple 
difference equations, they are not eigenfunctions of any difference operators. 
In addition, the corresponding variational problem is more difficult; while we 
can indeed solve the associated PDE, we have so far been unable to derive the 
solution. This is why we focus on a limiting case in the present paper. 
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Fix u±U2 — pqC, 2 , let ui, U2 = Q(y/p) as p — > 0, and similarly let C(i) ~ p 1 / 2 . 
In this limit, one has 

Km^(i,i) = C"(i) fc^-A 

or in other words the g-Racah weight(s) discussed above. In the notations of 
Sections 2-9, £ — K? - ^"* -1 " 2 . The corresponding limit of the elliptic MacMahon 
identity is 

g 2fc + l C 2 _ g t+J-l 
2-^1 11 „2fcA2 _ i+j 

ncaxbxc (i,j,fe)en 



^ _ gi+ j +fc -1^2 _ g i+j-fc-2) 

11 n _ a i+j+k-2\i(-2 _ i+j-k-i) ' 

l<i<a,l<j<b,l<k<c v ^ Mb ^ y 



or equivalently (performing the products over k and simplifying) 

V aim TT S-Zl —= TT 1^ 

y 11 (2 _ g i+]-c-2 11 1 - 

ncaxbxc l<i<a,l<i<b l<i<a,l<j"<6 y 



This, of course, becomes the usual MacMahon identity upon taking the limit 
£ — > oo, cf. Section 7.21 in [St]. 
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